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I . INTRODUCTION 

A. BACKGROUND 

The field of cyclic spectral analysis is growing in 
importance as non- stationary modulation schemes proliferate in 
modern communications. Traditionally, spectral analysis 
methods have assumed that signals of interest were stationary 
or had non- time -varying statistical parameters. This 
assumption is untrue for the cyclostationary signals created 
when using pulse or carrier amplitude modulation, phase or 
frequency carrier modulation or digital pulse or carrier 
modulation. Gardner [Ref. l:pp 419-453] discusses many of 
these signal types. Spread spectrum modulation techniques 
such as direct sequence PSK and frequency- hopped FSK also 
exhibit temporally varying spectral features [Ref. l:pp 453- 
457] . Reference 2 provides a tutorial on cyclic spectral 
analysis . 

Cyclic spectral analysis methods take advantage of 
spectral variation over time by correlating spectral 
estimations at discrete time intervals to locate signals that 
may not otherwise be identified. A direct sequence signal 
with a very wide bandwidth may be indistinguishable from a 
noisy background in the frequency spectral estimate of 
Figure 1 [Ref. 3:pp 2-9 - 2-10]. But this signal becomes 
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clearly identifiable in the cyclic spectral estimate of 
Figure 2 . 

One cyclic spectral technique is the frequency smoothed 
cyclic periodogram method (FSM) [Ref. l:pg 464] . This has 
been implemented in the Cyclic Spectral Analysis Software 
Package [Ref. 4] as the program sxaf. Hereafter, this 
software will be referred to as the SSPI package or the FSM. 
This method was used to generate Figures 1 and 2. 
Unfortunately, the FSM is less efficient than time smoothing 
methods for general spectral estimation [Ref. 5:pg 38]. 

B. THESIS GOALS 

The purpose of this thesis is to implement two time 
smoothing algorithms, the FFT Accumulation Method (FAM) 
[Ref. 5 : pp 44-47] and the Strip Spectral Correlation Algorithm 
(SSCA) [Ref. 5 : pp 47-48]. The FAM is less computationally 
complex than the FSM and the SSCA less so than the FAM. Both 
methods have been written to be compatible with the SSPI 
package. By integrating FAM and SSCA into the SSPI package 
future research may more easily be done in algorithm 
performance comparison and enhancement. 

Both the FAM and SSCA accept input data generated by the 
SSPI utility program pam. They also generate results 
compatible with the SSPI plotting utility plot_sxaf . Each 
program may optionally generate output in a format compatible 
with the MATLAB [Ref. 6] package which is in common use at the 
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Naval Postgraduate School. Input and output data are in ASCII 
file formats. 

Throughout this thesis the same BPSK signal data is used 
to illustrate the implementations of FAM and SSCA. A 
representative signal sample is plotted in Figure 3 . No noise 
was added to this signal. By estimating the frequency 
spectrum at different points in time, it is evident in Figure 
4 that the spectral features do indeed vary temporally. 

Chapter II describes the implementation of the FAM. 
Chapter III describes the implementation of the SSCA. Chapter 
IV describes a utility to estimate cycle frequencies for a 
given baseband frequency. The Sub-FFT Accumulation Method 
(SUBFAM) program [Ref. 7] gives a quick result for a specific 
subset of the spectral plane. Chapter V discusses specific 
algorithm performance comparisons. Finally, Chapter VI 
presents conclusions and recommendations for future areas of 
research. 
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o. 



Figure 



Figure 2 




1 Power Spectral Density of 3PSK Signal in Noise 
[Ref. 4:pg 18] 




in Noise [Ref. 4:pg 18] 
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Figure 3 Typical BPSK Signal Data 
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II. FFT ACCUMULATION METHOD 



A. THEORY 

The FFT Accumulation Method (FAM) [Ref. 5:pp 44-47] is a 
Fourier transform of correlation products between spectral 
components smoothed over time. Periodicities in the spectral 
components then become detectable. The cyclic spectral plane 
as shown in Figure 5 ranges in frequency from -f s /2 to +f s /2, 
and in cycle frequency from - f s to +f 5 where f s is the sampling 
frequency. For each unique (fj # a q ) cell in Figure 5, the FAM 
cross - spectral estimate is defined to be [Ref. 5:pg 44]: 



The FAM auto- spectral estimate is defined to be 
[Ref . 5 :pg 44] : 




-i2mq 

P 



( 1 ) 



c *o + ^« 

& XX T 



P-1 



( nL , fj) At =£ X t {iL, f k ) X? (rL, f t ) g c {n-r) e 



i = o 



-i2nrq 



(2) 



where: j is the average frequency bin (k+£)/2 

q is the difference or bandwidth frequency (k-£) 
g c (n-r) represents the Hamming window. 

The FAM was implemented by forming an array from x(kT) 
(0 < k < N-l) with rows which are N' points long from the 
input sample data. The starting point of each succeeding row 
is offset from the previous rows starting position by L 
samples. A Hamming window is applied across each row which is 
then Fast Fourier transformed and downconverted to baseband. 
The result at this point is a two-dimensional array with 
columns representing constant frequencies. Each column was 
point-wise multiplied in turn with the conjugate of the 
columns resulting from processing y(kT) (0 k <; N-l) . Each 
resultant product vector was Fast Fourier transformed and the 
low frequency half placed into the final cyclic spectral plane 
at the appropriate location in the cyclic spectral plane. 
Figure 6 shows a block diagram for the FAM cross- spectral 
estimate . 

The following subsections in this chapter discuss in 
detail what has been described in the previous paragraph. The 
cycle frequency resolution of FAM is Ace = f s /PL [Ref. 5:pg 44] 
where f s is the original data sampling rate. The frequency 
resolution of FAM is Af = tjN’ [Ref. 5:pg 44] . The time- 
frequency resolution product is AtAf = PL /N' [Ref. 5:pg 45] . 
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The command line format for calling the FAM program is 
provided in Appendix A. The FAM source code listing is 
provided in Appendix B. 

B . IMPLEMENTATION 

1. Input Channelization 

The input sample data is formed into a two-dimensional 
array. The array row length is equal to the number of input 
channels N ' . For a given number of input sample points N, a 
row size of N’ , and a chosen offset L, there are P = (N- 
N' ) / L ~ N/L rows formed. The choice of N' must take into 
consideration that ideally the time - frequency resolution 
product must be much greater than one [Ref. 5:pg 40]. N' 
should also be a power of 2 to avoid truncation or zero- 
padding in the FFT routines. L should be chosen to be less 
than or equal to N' /4 . 

The completely filled array is P rows by N' columns. 
Figure 7 shows how a small array is filled from a discretely 
sampled signal x(kT) when N' =16, P=8 and L=4 . The number 
inside each cell represents the value of k used to index on 
x(kT) to fill that location in the array. 

Figure 8 shows the magnitude of the original data for 
the example BPSK signal organized into the P by N' array where 
N=512 , N' =32, L=4 and P=128. The input data is assumed to be 
complex with a real and imaginary component to each sample 
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point. Figure 9 shows a single row of the array. The phase 
changes of the BPSK are evident . 

2 . Windowing 

A Hamming window [Ref 7:pg 467] is applied to each row 
of the array. The equation for the Hamming window is: 

w{r) =0.54-0.46cos( _^1LL ) , 0 <; r < N'-l (3) 

N -1 

A 32 -point Hamming window is plotted in Figure 10. It 
is applied to both the real and imaginary parts of the complex 
example array. The magnitude of the resultant array is shown 
in Figure 11. 

3. First FFT 

Each row of the windowed data array is Fast Fourier 
transformed to reveal the first spectral components. The 
resultant array is still indexed P rows by N' columns but now 
the column index relates to a specific bin of spectral 
frequencies. Figure 12 illustrates this relationship. 

Figure 13 shows the results of FFTing the BPSK example. 

4. Downconversion 

Each row of spectral components is downconverted to 
baseband through multiplication with the complex exponential, 

-i2nkmL 
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where: m is the row index, 0 s m < P-1 

k is the column index, 0 <; k <; N' - 1 

The magnitude of the exponential is unity over the 
array but the phase shows considerable variation. Figure 14 
shows the phase of the exponential over the P by N’ array. 
Figure 15 shows the phase for one representative row. The 
magnitude of the array remains unchanged from Figure 13. 

5. Multiplication 

For each cell in the cyclic spectral plane, one column 
of the array is multiplied with the conjugate of another 
column. The separation of the columns is determined by the 
desired frequency separation or cycle frequency (a=f k -f f ) . The 
mid-point between the columns is the frequency bin of interest 
(f = (f k +f f ) /2) . Figure 16 shows two columns to be conjugate 

multiplied for use in filling a specific f ,ot cell. Figure 17 
illustrates the f , o' bin values in the cyclic spectral plane 
for N‘ =4 . Figure 18 shows the corresponding two column 
indices used to form the product vector which is used to 
produce the cells of Figure 17. 

6 . Second FFT 

The product from the previous multiplication is FFT'd 
to yield a P-point result. Only the middle of the resulting 
spectrum is retained and stored into the designated t,ot cell. 
The upper and lower ends are undesireable because of increased 
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estimate variation at the channel pair ends [Ref. 5:pg 45]. 
Figure 19 illustrates which part of the estimate is retained 
and placed into the cell. Figure 20 shows the final result of 
the FAM auto- spectral process on the example BPSK signal. 

7 . Data Reduction 

The FAM program typically generates large output data 
files. For convenience, an option may be chosen to reduce the 
amount of output. By comparison sorting for the largest a 
value in an f , a cell, the number of ot slices output is reduced 
from N' P/ 2 to N' . Overall FAM program execution time 
increases accordingly to accomplish the sorts. 

Figure 21 illustrates the output full cyclic spectral 
plane storage array before sorting. Figure 22 shows the 
output array after data reduction has been completed. Figures 
23 and 24 plot the resulting spectral half -planes without and 
with data reduction respectively. 

Since all cells have an equal number of oi values, all 
sorts are of equal complexity. Further work on data reduction 
would require selection of a largest ot value from widely 
varying numbers of a values depending on the choice of N, N' 
and L. 

C . PERFORMANCE 

An established method of evaluating the complexity of an 
algorithm is to determine the number of floating point 
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operations that must be performed. For this estimate it is 
assumed that an auto- spectral estimate is being computed and 



the data is complex in every step. 


It is also assumed that 


the Hamming window coefficients 


and the downconversion 


coefficients have been previously 


computed and stored for 



later use. Each N-point FFT requires (N/2)*log 2 N complex 
floating point multiplies [Ref. 8:pg 506] or 4* (N/2 ) *log 2 N 
real floating point multiplies. The cost of any output data 
reduction is not considered here. 



Applying the window 


2*P *N' 


First FFT 


2*P*I\T *log,I\T 


Downconversion 


4*p*jV' 


Multiplication 


4 *N' *N' *P 


Second FFT 


2 *N' *N ' *P*log,P 



Total: (6+4*i\T ) *P*I\T + (2*P*2\T ) * (log 2 2\T + N' *log 2 P) (4) 

Since P-N/L (5) 

Total: (6+4*i\T ) *Ni\T /L + ( 2 *Ni\T / L ) ( 1 og 2 AT + N' log 2 N/L) (6) 
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Figure 5 Generic FAM Cyclic Spectral Plane 
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Figure 6 



Block Diagram of FAM Cross - spectral Estimate 
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Figure 7 Layout of A Sample Input Array Showing Input 
Sample Storage for N' = 16, N=48, L=4 and P=8. 
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Figure 9 Single Row of BPS K Signal Inpuc Array 
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Figure 11 



Example BPSK 



Signal Data After Windowing 



Co Tulin Index 



Scv *x 




H' -i 



frequency -f /2 
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t 1 2 -1 



Figure 12 Generic Array Row After First FFT 



13 




Figure 13 Example B?SK Signal Data After First FF 
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Figure 15 Phase For a Sinqle Row of Array 

Coluan Index 

0 1 - ► K' -1-15 * 




Figure 16 Generic Array Showing Two Columns to be 
Multiplied, f = -1, a = 2 



20 



i 




Figure 17 Generic Array Showing f , a values 




Figure 18 Generic Array Showing Column Index Values 
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Figure 19 Retained Section of the Second FFT Results 
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Figure 20 Final FAM Results of the BPSK Example, 
N=512 , AT* =32 and L=4 . 
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Figure 21 FAM Storage Before Data Reduction 



24 




Figure 22 FAM Storage After Data Reduction 
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Figure 23 FAM Result Without Data Reduction, 
N-512, N' =32 and L=4 . 
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Figure 24 FAM Result With Data Reduction, 
N-512, N' =32 and L=4 . 
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III. STRIP SPECTRAL CORRELATION ALGORITHM 



A. THEORY 

The Strip Spectral Correlation Algorithm (SSCA) 
[Ref. 5:pp 47-48] is a Fourier transform of correlation 
products between spectral and temporal components smoothed 
over time. Periodicities in the spectral components then 
become detectable. The cyclic spectral plane as shown in 
Figure 25 ranges in frequency from -f s /2 to +f s /2, and in cycle 
frequency from - f s to +f s where f s is the sampling frequency. 
For each unique (f j# a q ) strip in Figure 25, the SSCA cross - 
spectral estimate is defined to be [Ref. 5:pg 47] : 



, fjr+qAa 






qAa 



) Ac =£ X T^ r ' f k ) y* g(n-r) e 



-i2nrq 



(7) 



The SSCA auto- spectral estimate is defined to be 
[Ref . 5 : pg 47] : 



£ k +qAa 

Sxx T \ n » ' 
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where: k is a multiple of the frequency resolution, 

-N' / 2 s k s \N' / 2) -1 

q is a multiple of the cycle frequency resolution, 
- N' s q s (N‘ -1) 

g(n-r) represents the Hamming window. 

The SSCA was implemented by forming an array from x(kT) 
(0 < k < N-l) with rows which are N' points long from the 
input sample data. The starting point of each succeeding row 
is offset from the previous rows starting position by L 
samples. A Hamming window is applied across each row which is 
then Fast Fourier transformed and downconverted to baseband. 
The result at this point is a two-dimensional array with 
columns representing constant frequencies. Each row is 

transposed and replicated L times for a total of PL columns. 
Each column is then multiplied by a sample y(kT) (0 < k < (PL- 
1) ) . Each resultant row of the array is Fast Fourier 

transformed and placed into a strip of the final cyclic 
spectral plane at the appropriate location. Figure 26 shows 
a block diagram for the SSCA cross - spectral estimate. 

The following subsections in this chapter discuss in 

detail what has been described in the previous paragraph. The 
cycle frequency resolution of SSCA is A ot = f s /N [Ref. 5:pg 48] 
where f s is the original data sampling rate and N is the 

number of input samples used. The frequency resolution of 
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SSCA is equal to the sampling rate divided by the number of 
channels N' , A f = f % /N' [Ref. 5:pp 42,48] . The time- frequency 
resolution product is AtAf = N /N' [Ref. 5:pg 48] . 

The command line format for calling the SSCA program is 
provided in Appendix C. The SSCA source code listing is 
provided in Appendix D. 

B. IMPLEMENTATION 

1. Input Channelization 

The input sample data is formed into a two-dimensional 
array. The array row length is equal to the number of input 
channels N' . For a given number of input sample points N, a 
row size of N* , and a chosen offset L, there are P = (N- 
N‘ ) / L « N/L rows formed. The choice of N' must take into 
consideration that ideally the time- frequency resolution 
product must be much greater than one [Ref. 5:pg 40]. N' 
should also be a power of 2 to avoid truncation or zero- 
padding in the FFT routines. L should be chosen to be less 
than or equal to N' / 4. 

The completely filled array is P rows by N' columns. 
Figure 27 shows how a small array is filled from a discretely 
sampled signal x(kT) when N' = 16, P=8 and L=4. The number 
inside each cell represents the value of k used to index on 
x{kT) to fill that location in the array. 

Figure 28 shows the magnitude of the original data for 
the example BPSK signal organized into the P by N' array where 
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N=512, N' =32, L=4 and P=128. The input data is assumed to be 
complex with a real and imaginary component to each sample 
point. Figure 29 shows a single row of the array. The phase 
changes of the BPSK are evident. 

2 . Windowing 

A Hamming window [Ref 7:pg 467] is applied to each row 
of the array. The equation for the Hamming window is: 

w(r) =0.54-0.46cos( JiILil ) t Q z r <, N'-l (9) 

N'-l 

A 32 -point Hamming window is plotted in Figure 30. It 
is applied to both the real and imaginary parts of the complex 
example array. The magnitude of the resultant array is shown 
in Figure 31. 

3 . First FFT 

Each row of the windowed data array is Fast Fourier 
transformed to reveal the first spectral components. The 
resultant array is still indexed P rows by N* columns but now 
the column index relates to a specific bin of spectral 
frequencies. Figure 32 illustrates this relationship. 

Figure 33 shows the results of FFTing the BPSK example. 

4. Downconversion 

Each row of spectral components is downconverted to 
baseband through multiplication with the complex exponential, 
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where: m is the row index, 0 ^ m s P-1 

k is the column index, 0 s k s 2\T - 1 

The magnitude of the exponential is unity over the 
array but the phase shows considerable variation. Figure 34 
shows the phase of the exponential over the P by N' array. 
Figure 35 shows the phase for one representative row. The 
magnitude of the array remains unchanged from Figure 33. 

5. Replication 

Each row is copied into one column of an empty N' by 
PL array. It is then replicated in the L-l adjacent columns. 
Figure 36 illustrates this process. 

6. Multiplication 

Each column of the array is pointwise multiplied with 
the conjugate of a sample value y(kT) . There are PL«N columns 
and PL samples from y(kT). Figure 37 illustrates the 
conjugate multiplication process. 

7 . Second FFT 

Each row from the previous multiplication is Fast 
Fourier transformed to yield a PL-point result. Each 
resulting vector is stored into a strip of the cyclic spectral 
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plane. Figure 3 8 shows the final result of the SSCA auto- 
spectral process on the example BPSK signal. 

8. Data Reduction 

The SSCA Program typically generates large output data 
files. For convenience, an option may be chosen to reduce the 
amount of output. By comparison sorting for the largest a 
value in an f,a cell, the number of a slices is reduced from 
N to N/L. Overall SSCA execution time increases accordingly 
to accomplish the searches. 

Figure 39 illustrates the output full cyclic spectral 
plane storage array before sorting. Figure 40 shows the 
output array after data reduction has been completed. Figures 
41 and 42 plot the resulting spectral half -planes without and 
with data reduction respectively. 

Since all cells have an equal number of ot values, all 
sorts are of equal complexity. Further work on data reduction 
would require selection of a largest a value from widely 
varying numbers of a values depending on the choice of N, N' 
and L. 

C . PERFORMANCE 

An established method of evaluating the complexity of an 
algorithm is to determine the number of floating point 
operations that must be performed. For this estimate it is 
assumed that an auto- spectral estimate is being computed and 
the data is complex in every step. It is also assumed that 
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the Hamming window coefficients and the downconversion 
coefficients have been previously computed and stored for 
later use. Each N-point FFT requires (N/2)*log : N complex 
floating point multiplies [Ref. 8:pg 506] or 2*N*log 2 N real 
floating point multiplies. The cost of any output data 
reduction is not considered here. 



Applying the window 2*P*i\T 

First FFT 2*P*N' *log 2 tf' 

Downconversion 4*P*N' 

Multiplication 4*N *N' 

Second FFT 2*N*N' *log 2 N 



Total : 


2*N' * ( 6*P + 4 *N + (2*P + 2 *N) 


*log 2 N) 


(10) 


since 


P«N/L 




(ID 


Total : 


2*i\T * ( (6*N/L + 4 *N) + (2*N/L 


+ 2*N) *log 2 N) 


(12) 
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Figure 25 Generic SSCA Cyclic Spectral Plane 
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Figure 27 Layout of A Sample Input Array Showing Inpuu 
Sample Storage for N' =16, N=48, L=4 and P-8. 
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Figure 28 Example BPSK Signal Data Array, N' = 32, P=123, L=4 
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Figure 31 Example BPSK Signal Data After Windowing 
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Figure 32 Generic Array Row After First FFT 
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Figure 33 Example 3PSK Signal Data After First FFT 




Figure 34 Phase of the Downconversion Exponential 
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Figure 35 Phase For a Single Row of Array 
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Figure 36 Replication Process 
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Figure 37 Multiplication of Columns by y* 
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Figure 38 Final SSCA Results of the BPSK Example, 
N=512, N' =32 and L=4 . 
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Figure 39 SSCA Storage Array Before Data Reduction 



46 




Figure 40 SSCA Storage Array After Data Reduction 
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Figure 41 BPSK Example Results Without Data Reduction, 
N=5 12 , N' - 32 and L=4 . 
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Figure 42 BPSK Example Results With Data Reduction, 
N=512 , N' =32 and L=4 . 
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IV. SUB-FFT ACCUMULATION METHOD 



A. THEORY 

In the course of this thesis work a need arose for a 
flexible and customized Spectral Correlation Function (SCF) 
[Ref. 7] to determine the presence of phase- shift keyed 
signals in experimental data. The SCF [Ref. 7] is 

represented : 

H-l r - ; 

S£y(f) = £ x(r)yir) e OSF x OSF y e N (13) 

r = 0 



where : 

N is the number of data samples used 

f, is the intermediate modulation frequency 

f s is the sampling frequency 

OSF is a one-sided bandpass filter 

The program written to implement this algorithm is named 
the Sub-FFT Accumulation Method (SUBFAM) program to 
distinguish it from the more general purpose FAM program. The 
command line format for calling the SUBFAM program is provided 
in Appendix E. The SUBFAM source code is provided in 
Appendix F. 
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B. IMPLEMENTATION 



Figure 43 illustrates the operations performed in the 
SUBFAM program. A set of N points is obtained from within a 
file of signal data values. Each set of N points is processed 
through a one-sided bandpass filter to remove either all 
positive or negative frequencies as desired. Downconversion 
of each file to baseband from the original sampling and 
transmission frequency bands follows. The results are then 
correlated through a complex conjugate multiplication. A 
final N-point Fast Fourier Transform yields the spectral 
correlation result. 

C . PERFORMANCE 

The SUBFAM program performed as required. Figure 44 
illustrates the results of processing test data containing 
phase- shift -keyed signal samples. The peaks at the chip 
frequencies are apparent. Figure 45 shows the results of 
correlating test data with data containing only white noise. 
The lack of correlation between spectral features yields 
results which are four orders of magnitude less than in the 
successful signal detection of Figure 44. 
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Figure 43 SUBFAM Program Block Diagram 
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Figure 44 SUBFAM Program Result With Signal, N=4096 
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Figure 45 SUBFAM Program Result With Noise, N=4096 



54 



V. ALGORITHM PERFORMANCE 



The complexity of the FAM and SSCA algorithms is estimated 
in Chapter II Section C and Chapter III Section C 
respectively. Both require on order of Nlog : N floating point 
multiplies where N is the original number of sample points 
used. Brown and Loomis [Ref. 9] explore the complexity of the 
FAM, SSCA and FSM algorithms in some detail. Their results 
are consistent with an order Nlog : N multiplies for FAM and 
SSCA. 

The cyclic spectral plane is symmetric about the f=0 and 
the ce=0 axes for real signals. For the most efficient 
performance it is only necessary to compute one quadrant of 
the cyclic spectral plane. To compare the performance of FAM 
and SSCA with FSM, complexity figures derived from Reference 
9 for quadrant complexity will be used. The three algorithms 
require the following numbers of floating point multiplies: 

FAM Complexity [Ref. 9:pg 18] 

2Ni\7' / log 2 ( ) + 8Nlog 2 (N / ) + 4NiV / = 20N (14) 

N 
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SSCA Complexity [Ref. 9:pg 20] 



NiV / log 2 (N) + 2N N' + 8Nlog Z (N / ) + 12N (15) 



FSM Complexity [Ref. 9:pg 16] 

N 2 + 2Nlog 2 (N) (16) 

The complexity of four example cases with varying N and N ' 
are estimated below using equation (14) , (15) and (16) . 



N 


N' 


FAM 


SSCA 


FSM 


1024 


64 


5 . 3 * 10 5 


3 . 6 * 10 5 


10 6 


2048 


128 


2 . 0*10 6 


1 .4*10 6 


4 . 2 * 10 6 


32,768 


512 


1 . 5 *10 8 


1.1*10 8 


10 9 


1, 048,576 


8192 


8 . 0*10 10 


7. 0*10 10 


10 12 



It is evident that FAM and SSCA require substantially 
fewer multiplies than FSM. Particularly with larger values of 
N , FAM and SSCA require Nlog 2 N while FSM requires N 2 floating 
point multiplies. 
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VI. CONCLUSIONS 



A. SUMMARY 

The purpose of this thesis was to implement the FAM and 
SSCA for use in cyclic spectral analysis research. Both 
algorithms were successfully implemented in a manner 
compatible with the SSPI analysis package to facilitate future 
work. By allowing the user a choice of two output file 
formats, results may easily be used in either MATLAB [Ref. 6] 
or SSPI [Ref. 4] for further signal processing applications. 

It has been shown that the FAM and SSCA algorithms 
generate results consistent with FSM but require substantially 
fewer floating point multiplies than FSM. This is especially 
true as the number of sample data points used increases. 

B. RECOMMENDATIONS 

In the course of this thesis it has become apparent that 
the FAM and SSCA algorithms consist of inherently parallel 
operations. Roberts and Loomis explore these possibilities in 
Reference 10. The high-level sequential nature of these 
algorithms also lend them to pipelining architectures. By 
combining a parallel or multi -processor approach with 
pipelining a large improvement in performance should be 
obtainable. This would be espectially true in applications 
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requiring large sets of signal data samples or repeated, 
time- sequenced cyclic spectral plane estimations. 

Two areas which could benefit from further work are 
automatic signal identification and signal parameterization 
algorithms. The results from cyclic spectral analysis 
algorithms such as FAM and SSCA may be input into Linear 
Associator or Mapping neurocomputer networks as described in 
Reference 11. Intensity transformations using digital image 
processing techniques from Reference 12 also appear to have 
promise. Both classes of techniques could have utility in 
automatic signal identification and parameterization. They 
also have the added benefit that they lend themselves to 
parallel and pipelined computing architectures. 
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APPENDIX A. FAM PROGRAM USE 



Correct commands are: 
for the cross -f am 

fam inputfilel inputfile2 outputfile n nprime 1 Oflag 
fam inputfilel inputfile2 outputfile n nprime 1 Oflag red 

or for the auto- fam 

fam inputfilel outputfile n nprime 1 Oflag 

fam inputfilel outputfile n nprime 1 Oflag red 

where: inputfilel is the file containing the x signal samples 

inputfile2 is the file containing the y signal samples 

outputfile is where the spectrum values will be placed 
n is the number of samples to use from the inputfiles 
and it must be a power of 2 
nprime is the group size of the input datasets and it 
must be a power of 2 

1 is the starting offset of subsequent datasets and 
it must be a power of 2 

Oflag indicates the output filetype, ascii or binary 
-asc for an ascii, MATLAB compatible file, 
full cyclic spectral plane 
-plo for a plot_sxaf , SSPI compatible file, 
half cyclic spectral plane 
red reduces the amount of data output 

note: the first two entries in every input file are expected 

to be the datatype and n. datatype = 1 for real 
values, 2 for complex values. n is the number of 
samples contained in the file, one sample per line. 

input file format: 

datatype n 
sample #1 
sample #2 



sample #n 
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output file ASCII MATLAB format: 



nprime n 

output (1) (1) 
output (1) (2 ) 



output (nprime ) (n) 



output file SSPI plot_sxaf format: 

datatype (1 for real, 2 for complex) 
number_of_alphasalpha_minalpha_max 
number_o f _f reqs f req_minf req_max 

value of alpha #1 

number of freqs in #lf req_minf req_max 
spectrum at freq_min 



spectrum at freq_max 



value of alpha #alpha_max 

number of freqs in #alpha_maxf req_minf req_max 
spectrum at freq_min 



spectrum at freq_max 
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APPENDIX B. FAM PROGRAM LISTING 



61 



Dec 9 12:21 1992 fam/fam.c Page 1 



^include <stdlib.h> 

^include <stdio.h> 
i include <math.h> 

$ include " /home3/carter/thesis/SSPI/pam/ f ft . c" 

# include " /home3/carter/thesis/SSPI/pam/radix . c" 
main (argc, argv) 
int argc; 
char *argv[]; 

{ 

COMPLEX *x, **s3, **y3, **dwnconv, **s4; 

COMPLEX *tempfft, *yl, *sl, **s2, **y2, ‘window; 

/* x passes data to/from fft routine 
s3 holds results of FFTing s2 
y3 holds results of FFTing y2 
dwnconv holds downconversion coefficients 
s4 holds correlation multiply results 
tempfft passes data to/ from fft routine 
yl holds y samples from input file 2 
si holds x samples from inputfile 1 
s2 holds channelized x samples 
y2 holds channelized y samples 
window holds Hamming window values 

*/ 

int i, j, k, type, n, file_n, nprime, p, 1, direction, norm; 
int a, b, c, cross, num_f , max_num_f , da t a_t yp e , max_num_a 1 f ; 
int tempi, temp j , halfp, reduce, curr_alf ; 
int * out int; 

float numreal, numimag, convfac, mainfac, num_alf, fjmin, f_max; 

float alpha_min, alpha_max, f_min_all, f_max_all; 

float tempreal, tempimag, bigmag, tempmag; 

float twopi=6. 28318530718; 

float *outreal; 

double z,y; 

char *infilel, + infile2, *outfile; 

FILE + ifpl, *ifp2, *ofp; 

/* check for the correct number of input arguments 
the correct commands are: 

fam inputfilel inputfile2 outputfile n nprime 1 Oflag reduce 
or 

fam inputfilel inputfile2 outputfile n nprime 1 Oflag 
or 

fam inputfilel outputfile n nprime 1 Oflag reduce 
or 

fam inputfilel outputfile n nprime 1 Oflag 

where: inputfilel is where the x signal samples are expected to be 

inputfile2 is where the y signal samples are expected to be 

outputfile is where the spectrum values will be put 

n is the number of samples to use from the inputfiles 
nprime is the group size of input datasets 
1 is the starting offset of subsequent datasets 
Oflag indicates the output filetype, ascii or binary 

-asc for an ascii, matlab compatible file, full plane 
-plo for a plot_sxaf compatible file, half plane 
reduce is an option to reduce the amount of output 
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note: datatype and n, are expected to ce at the top of the input file 

datatype = 1 for real, 2 for complex 
n is the number of samples following 

written by LCDR Nancy J. Carter, USN NPGS Monterey CA 01 October 1992 
to implement the Frequency Accumulation Method of cyclic spectral analysis 
as described in the Brown/Gardner/Loomis IEEE paper 

20 October 1992 .... changed to ascii MATLAB output format 

22 October 1992 .... added plot_sxaf output format compatible with SSPI 

23 November 1992 ... .added option to reduce amount of output 

* / 

CHECK INPUT VALUES */ 



/* */ 
if ( (argc ! = 7)&&(argc ! = 8)&&(argc 1= 9) ) { 
pnntf ( " fam. . . . fatal errorin" ) ; 

print f(" incorrect number of calling arguments \n" ) ; 

print f(" correct formats are : in"); 

printfC fam inputfilel mputfile2 output file n nprime 1 Oflag reducein") 

pnntf (" fam inputfilel mputfile2 outputfile n nprime i Oflagin"); 

printf(" fam inputfile outputfile n nprime 1 Oflag reducein"); 

printf(" fam inputfile outputfile n nprime 1 Oflagin"); 

printfC wherein"); 

printfC inputfilel contains the x samplesin"); 

printfC inputfile2 contains the y samplesin"); 

printfC' outputfile will contain the resultsin"); 

printfC' n is the number of input samples to use\n") ; 

printfC' nprime is the group size of input datasetsin" ) ; 

printfC' 1 is the offset of subsequent datasetsin"); 

printfC' Oflag is the output file formatin''); 

printfC' Oflag = -asc, an ascii file is producedin" ) ; 

printfC' Oflag = -plo, a plot_sxaf file is producedin"); 

printfC' reduce is an option for reduced amount of output in" ) ; 

exit () ; 

} 

if (argc==7) { /* this is an auto- fam, no output reduction */ 

crossed ; 



inf ilel=argv [ 1 ] ; 
outfile=argv [2] ; 
n=atoi (argv [3] ) ; 
nprime=atoi (argv [4 ] ) ; 
l=atoi (argv [5] ) ; 
reduce=0; 

} 

if ( (argc==8) && (argv [7 ] [0] 1 =' r ' ) ) { /* this is a cross-fam, no output reduction */ 

cross=l; 

inf ilel=argv [ 1 ] ; 
inf ile2=argv [2 ] ; 
out f i le=ar gv [ 3 ] ; 
n=atoi (argv [4] ) ; 
nprime=atoi (argv [5] ) ; 
l=atoi (argv [ 6] ) ; 
reduce=0; 

} 

if ( (argc==8) && (argv [7 ] [0] ==' r' ) ) { /* this is an auto- fam, with output reduction */ 

cross=0; 
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mfilel=argv [1] ; 
outf ile=argv [2 ] ; 
n=atoi (argv [ 3 ] ) ; 
npnme=atoi (argv [4 J ) ; 
l=atoi (argv [5] ) ; 
reduce=l; 

} 

if (argc==9) { /* this is a cross-fam, with output reduction */ 

cross=l ; 

inf ilel=argv [ 1 ] ; 
inf ile2=argv [2] ; 
outfile=argv [3] ; 
n=atoi (argv [ 4 ] ) ; 
npnme=atoi (argv [ 5 ] ) ; 
l=atoi (argv [ 6] ) ; 
reduce=l; 

} 

/* verify that n, nprime and 1 are powers of 2 */ 

i=n%2; 

if (l 1 =0) { 

print f ( " f am .... fatal error\n " ) ; 

pnntf (" calling argument n is not a power of 2\n"); 

exit () ; 

} 

i=nprime%2; 



pnntf (" fam. . . . fatal error\n" ) ; 

printf (" calling argument nprime is not a power of 2\n"); 



exit () ; 

} 

j=l%2; 
if (j!=0) { 



printf ("fam. . . . fatal error\n" ) ; 

printf (" calling argument 1 is not a power of 2\n"); 

exit 0 ; 

} 

/* */ 

/* open input file 1 , prepare to get the x signal sample data */ 
/* V 



ifpl = f open (inf ilel, "r" ) ; 

fscanf (ifpl, "%i %i", &type, &file_n) ; 

/* verify that file_n is greater than or equal to n */ 
if (file_n<n) { 



printf ("fam. . . . fatal error\n" ) ; 

printf (" inputfilel does not contain enough samples\n"); 

fclose (ifpl) ; 
exit ( ) ; 

} 

/* find p - the number of datasets (nprime long) */ 
p= ( (n-nprime) /l) ; 

/* */ 
/* */ 



/* READ IN DATA SAMPLES */ 

/* */ 

/* allocate space to read in 
/* */ 



the sample values */ 
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si* (COMPLEX*) cailoc (n, sizeof (COMPLEX) ) ; 
if (sl==NULL) { 

pnntf (" fam. ... insufficient space to allocate si\n M ); 
exit () ; 



if (type==l) { 

/ * read in the real sample values */ 

for (i=0; i < n; i++) { 

fscanf (ifpl, "%e\n", snumreal) ; 
sl[i] . r=numreal; 
si [i] . i=0 . 0; 

} 

else { 

/* 

/* read in the complex sample values */ 

/* 

for ( 1=0 ; i < n; i++) { 

fscanf (ifpl, "%e %e\n", snumreal, Snumimag) ; 

sl[i] .r=numreal; 
si [i] . i=numimag; 

} 

} 

/* */ 

/* close the input file #1 */ 

/* */ 

fclose (ifpl) ; 

/* */ 

/* if this is a cross-fam go get y samples */ 

/* */ 

if (cross==l) { 

/* */ 

/* open input file2, prepare to get the y signal sample data */ 
/* */ 

ifp2 = fopen (infile2, "r" ) ; 
fscanf (ifp2 , "%i %i", &type, &file_n) ; 

/* verify that file_n is greater than or equal to n */ 
if (file_n<n) { 

printf ( " fam. . . . fatal errorXn" ) ; 

pnntf (" input filel does not contain enough samples\n") 

fclose (ifpl) ; 
exit ( ) ; 

} 

/* */ 

/* READ IN Y DATA SAMPLES */ 

/* */ 

/* allocate space to hold the sample values */ 

/* */ 

yl* (COMPLEX* )calloc(n, sizeof (COMPLEX) ) ; 
if (y 1==NULL) { 

print f (” fam. ... insufficient space to allocate yl\n"); 
exit () ; 

} 

*/ 



b5 



/ 
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* read in the complex sample values * 

for ( 1=0 ; i < n; i++) { 

fscanf (ifp2, " %e %e\n" , &numreal, inumimag) ; 
y 1 [ i ] . r=numreal; 
yl[i] .i=numimag; 

} 

/* */ 

* close the inputfile2 */ 

*/ 

fclose (ifp2 ) ; 

} /* end if cross=l */ 

/* */ 

/* FORM DATASETS */ 

/* */ 

/* allocate space to hold the p-by-nprime datasets */ 

/* */ 

s2= (COMPLEX**) calloc (p, sizeof (COMPLEX*) ) ; 
if (s2==NULL) { 

printf ("fam. ... insufficient space to allocate s2\n"); 
exit ( ) ; 

} 

for (i=0; i<p; i++) { 

s2 [i ] = (COMPLEX* ) calloc (npnme, sizeof (COMPLEX) ) ; 
if (s2 [i ] ==NULL) ( 

printf ("fam. ... insufficient space to allocate s2\n"); 
exit ( ) ; 

} 

} 

/* V 

/* form p datasets of nprime samples from the original sample stream 
/* */ 

for (i=0; i<p; i++) { 

for (j=0; j<nprime; j++) { 
k=i*l+j; 
s2[i] [ j ] =sl [k] ; 

} 

} 

/* */ 

/* if this is a cross-fam form y datasets */ 

/* */ 

if (cross==l) { 

/* */ 

/* */ 

/* allocate space to hold the p-by-npnme datasets */ 

/* */ 

y2= (COMPLEX* * ) calloc (p, sizeof (COMPLEX* ) ) ; 
if (y2==NULL) { 

printf ("fam. ... insufficient space to allocate y2\n"); 
exit ( ) ; 

} 

for (i=0; i<p; i++) { 

y2 [ i ] = (COMPLEX* ) calloc (nprime, sizeof (COMPLEX) ) ; 
if (y2 [i] —NULL) { 

printf ("fam. ... insufficient space to allocate y2\n"); 
exit ( ) ; 
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} 

* * / 

* form p datasets of nprime samples from the original sample stream * 

for (i=0; i<p; i++) { 

for (j=0; j<nprime; j++) { 
k=±*l+j; 
y2[i] [j]=yl[k]; 

} 

} 

} /* end if cross=l */ 

/* */ 

/* APPLY WINDOW TO DATASETS */ 

/* */ 

/* allocate space to hold the window multiplicand values */ 

/* */ 

window= (COMPLEX* ) calloc (nprime, sizeof (COMPLEX) ) ; 
if ( window==NULL ) { 

printf ("fam. ... insufficient space to allocate sl\n"); 
exit ( ) ; 

} 

/* */ 

/* calculate the window values for the "nprime'' sample wide window */ 
/* */ 

for (i=0; i<nprime; i++) { 
y= (twopi*i) / (nprime-1) ; 
z=cos (y) ; 

window [i] .r= 0.54 - (0.46*z); 
window [i] . i= window [i] .r; 

} 

/* */ 

/* apply the window to all rows of s2 */ 

/* */ 

for (i=0; i<p; i++) { 

for (]=0; j<npnme; j++) { 

s2 [i] [ j] .r=s2 [i] [ j] . r*window [ j ] .r; 
s2 [i] [ j] .i=s2 [i] [ j] .i*window[^] . i; 

} 

} 

/* */ 

/* if this is a cross-fam apply window to y2 */ 

/* */ 

if (cross==l) { 

/* */ 

for (i=0; i<p; i++) { 

for (j=0; j<nprime; j++) { 

y2 [i] [ j] • r=y2 [i] [ j] .r*window[ j] .r; 
y2 [i] [ j] . i=y2 [i] [ j] . i*window[ j] .i; 

} 

} 

} /* end if cross=l */ 

/* */ 

/* FFT EACH DATASET ROW */ 

/* */ 

/* allocate space to hold rows of data for passing to the FFT routine 
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x= (COMPLEX* ) calloc (npnme, sizeof (COMPLEX) ) ; 
if (x==NULL) { 

printf (" fam. ... insufficient space to allocate x\n"); 
exit () ; 

} 

*/ 

* allocate space to hold rows of results from the FFT routine */ 
, * */ 

s3= (COMPLEX**) calloc (p, sizeof (COMPLEX*) ) ; 
if (s3==NULL) { 

printf ("fam. ... insufficient space to allocate s3\n" ); 
exit () ; 

} 

for (i=0; i<p; i++) { 

S3 [ l ] = (COMPLEX* ) calloc (npnme, sizeof (COMPLEX) ) ; 
if (s3 [i] ==NULL) { 

printf ("fam. ... insufficient space to allocate s3\n"); 
exit ( ) ; 

} 

} 

/* */ 

/* take an FFT of each row of s2 */ 

/* */ 

/* set values in arguments sent to fft */ 

directional; 

norm=l; 

for (i=0; i<p; i++) { 

/* copy row of s2 into complex array */ 

for (j=0; j<nprime; j++) { 
x [ j ] =s2 [ i ] [ j ] ; 

} 

/* go get FFT performed */ 

fft (x, nprime, direction, noon) ; 

/* copy x values into a row of s3 */ 

for ( j— 0 ; j<npnme; j++) { 
s3 [i] [ j]=x[ j] ; 

} 

} 

/* */ 

/* if this is a cross-fam take an FFT of each row of y2 */ 

/* */ 

if (cross==l ) { 

/* */ 

/* */ 

/* allocate space to hold rows of results from the FFT routine */ 
/* */ 

y3= (COMPLEX**) calloc (p, sizeof (COMPLEX*) ) ; 
if (y 3==NULL ) { 

printf ("fam. ... insufficient space to allocate y3\n"); 
exit ( ) ; 

} 

for (i=0; i<p; i++) { 

y3 [i] = (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 
if (y3 [i] ==NULL) { 

printf ("fam. ... insufficient space to allocate y3\n"); 
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exit ( ) ; 



} 

for (i= 0 ; i<p; i++) { 

/* copy row of y 2 into complex array */ 

for (j= 0 ; jcnprime; 3 ++) { 
x[j]=y 2 [i] [ j] ; 

} 

/ * go get FFT performed * 

f ft (x, npnme, direction, norm) ; 

/* copy x values into a row of y3 */ 

for (j= 0 ; ]<nprime; j++) { 
y 3 [ i ] [ j ] =x [ j ] ; 

} 

} 

} /* end if cross=l */ 

/* */ 

/* */ 

/* DOWNCONVERT EACH ROW */ 

/* */ 

/* allocate space to hold the downconversion multiplicands */ 

/* */ 

dwnconv= (COMPLEX^*) calloc (p, sizeof (COMPLEX*) ) ; 
if (dwnconv==NULL) { 

printf (" fam. ... insufficient space to allocate dwnconv\n"); 
exit () ; 

} 

for (i= 0 ; i<p; i++) { 

dwnconv [ i] = (COMPLEX* ) calloc (nprime, sizeof (COMPLEX) ) ; 
if (dwnconv [i] —NULL) { 

printf ("fam. .. .insufficient space to allocate dwnconv\n"); 
exit () ; 

} 

} 

/* */ 

/* downconvert each of the transform sequences */ 

/* */ 

/* calculate the down conversion multipliers */ 

/* */ 

mainfac=twopi*l/nprime; 

for (i= 0 ; i<p; 1 ++) { 

for ( j= 0 ; jcnprime; j++) { 
convfac=i* j*mainfac; 
dwnconv [ 1 ] [ j ] . r=cos (convfac) ; 
dwnconv [ i] [ j ] . i= (- 1 . 0 ) *sin (convfac) ; 

} 

} 

/* multiply downconversion factor against each frequency value */ 
for (i= 0 ; i<p; i++) { 
for ( j* 0 ; j<nprime; j++) { 

numreal=s3 [i] [j] . r*dwnconv [i] [j] .r - s3[i] [j] . i*dwnconv [i] [j] . 1 ; 
s3 [i] [ j ] . i=s3 [i] [ j ] . r*dwnconv [i] [ j] .i + s3 [i] [ j] . i*dwnconv [i] [ j ] . r 
s3 [i ] [ j ] . r=numreal; 

} 

} 

if (cross==l ) { 
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* */ 

* if this is a cross-fam downconvert y3 */ 

* multiply downconversion factor against each frequency value */ 
for (i= 0 ; i<p; i++) { 

for (j= 0 ; j<nprime; j++) { 

numreal=y3 [i] [ 3 ] . r*dwnconv [i] [ j ] . r - y3[i] [j] . i*dwnconv [i] [j] . i; 
y3 [i] [ j] . i=y3 [i] [ j] . r*dwnconv [i] [ j ] . i + y 3 [ i ] [ j ] . i*dwnconv [i ] [ j] . r; 
y 3 [ i ] [j] . r=numreal; 

} 

} 

/* end if cross=l */ 

*/ 

/* MULTIPLY COLUMNS */ 

/* */ 

/* allocate space to hold the correlation multiply (column multiply) results */ 

/* */ 

a=nprime*p/ 2 ; 

s4= (COMPLEX* * ) calloc (a, sizeof (COMPLEX* ) ) ; 
if (s4==NULL) { 

printf ("fam. ... insufficient space to allocate s4\n"); 
exit () ; 

} 

for (i= 0 ; i<a; i++) { 

s4 [i]= (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 
if (s4 [i] ==NULL) { 

printf (" fam. ... insufficient space to allocate s4\n"); 
exit ( ) ; 

} 

} 

/* */ 

/* allocate space to hold columns of s4 for passing to the FFT routine */ 

/* */ 

tempf ft= (COMPLEX*) calloc (p, sizeof (COMPLEX) ) ; 
if (tempf ft ==NULL) { 

print f (" fam. ... insufficient space to allocate tempf ft %i by l\n",p); 
exit () ; 

} 

/* set values in arguments sent to fft */ 

direct ion=l ; 

norm=l; 

if (cross== 0 ) { /* auto-fam */ 

/* */ 

/* multiply columns of s3 */ 

/* */ 

for (a=nprime-l; a> - 1 ; a — ) { 
for (b= 0 ; b<nprime; b++) { 

/* multiply the columns of s3 with other columns of s3 */ 
for (j= 0 ; j<p; j++) { 

tempf ft [ j ] . r= (s3 [ j ] [a] .r*s3[j] [b] .r) + (s3[j] [a] .i*s3[j] [b] . 1 ); 
tempf ft [ j ] . i= (s3 [ j ] [a] . i*s3 [ j ] [b] . r) - (s3 [ j ] [a] .r*s3[j] [b] . 1 ); 

} 

/* */ 

/* FFT EACH COLUMN */ 

/* */ 

/* go get FFT performed */ 
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f ft (tempf ft , p, direction, norm) ; 

/ * copy tempf ft result values back into a row of s4 
for ( j= (p/4) ; j< ( (3*p/4) -1) ; j++) { 
c= ( (npnme-a-1) *p/2) + 3 - (p/4) ; 
s4[c] [b] =tempf ft [ j ] ; 

} 

} /* end for b=0 . . . */ 

} /* end for a=0 ... * / 

} /* end auto-fam */ 

else { /* cross-fam */ 

/* */ 

/* multiply columns of s3 with y3 V 
/* */ 

for (a=nprime-l; a> -1; a — ) { 

for (b=0; b<npnme; b++) { 

/* multiply the columns of s3 with other columns of y3 
for ( j=0; 3 <p; ]++) { 

tempf ft [ 3 ] . r= (s3 [ j ] [a] .r*y3[;j] [b] .r) + (s3[;j] [a] .i*y3[;j] [b] . 1 ) ; 
tempf ft [j] . 1 = (s3 [ j ] [a] . i*y3 ( j] [b] . r ) - (s3 [ j ] [a] . r*y3 [ j] [b] . 1 ) ; 

} 

/* */ 

/* FFT EACH COLUMN */ 

/* */ 

/* go get FFT performed * / 

f ft (tempf ft, p, direction, norm) ; 

/* copy tempfft result values back into a row of s4 */ 

for (j=(p/4); j< ( (3*p/4) -1) ; j++) { 
c— ( (nprime-a-1) *p/2) + j- (p/4) ; 
s4 [c] [b] =tempf ft [ j ] ; 

} 

} /* end for b=0 ... */ 

} /* end for a=0 ... */ 

} /* end cross-fam */ 

/* */ 

/* OUTPUT */ 

/* */ 

halfp=p/2; 

if (reauce==0) { /* do not reduce the amount of output */ 

if ( ( (argc==8)&& (argv[7] [ 1 ] ==' a' ) ) | | ( (argc==7)&& (argv[6] [l]= =, a' ) ) ) { 

/* make an ascii output file of the whole spectral plane */ 

/* open the output file and place header information in it */ 
ofp = fopen(outfile / "w"); 
max_num_alf=nprime*p/2 ; 
max_num_ f = np r ime ; 

fprintf (ofp, " %i %i\n" , max_num_alf , max_num_f ) ; 
for (i=0; i<max_num_alf ; i++) { 
for (j=0; j<nprime; j++) { 

fprintf (ofp, "%e %e\n", s4 [i] [j] . r, s4 [i] [ j] . i) ; 

} /* end for j=. . . */ 

} /* end for i=. . . */ 

fclose (ofp) ; 

} /* end if ascii... */ 

else { /* make a plot_sxaf no reduced output file */ 
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/* open the output file and place header information in it */ 
ofp = fopen (outf ile, "w" ) ; 

/* place header information in the file */ 
data_type= 2 ; /* real_or_complex */ 
max_num_alf=nprime*p/ 2 ; /* num_alpha */ 

alpha_min= 0 . 0 ; /* alpha_min */ 

alpha_max=l . 0 ; /* alpha_max */ 

max_num_f =npnme ; /* max_num_f */ 

f_min_all= -0.5; / * f_min_all '/ 

f_max_all= 0.5; /* f_max_all */ 

fprintf (ofp, "%i %i %e %e %i %e %e\n" , data_type, max_num_alf , alpha_min, 
alpha_max, max_num_f , f_min_all, f_max_all) ; 
for (i= 0 ; icnpnme; i++) { /* group of lines */ 

for (j= 0 ; j<(p/ 2 ); j++) { /* line in the group */ 

/* place alpha line header information in file */ 
curr_alf= (i*p/ 2 ) + 3 ; /* alpha number */ 

num_f=npnme-i; /* num_f */ 

f_min= - 0 . 5+ ( 0 . 5*i/nprime) ; / * f_min */ 
f_max= 0 . 5- ( 0 . 5*i/nprime) ; /* f_max */ 
fprintf (ofp, "%i %i %e %e\n" , curr_alf , num_f , f_min, f_max) ; 
for (k=i; k<nprime; k++) { /* points on the line */ 

tempi= (p/ 2 * (nprime- 1 ) ) - ( (k-i) *p/ 2 ) +j ; 
temp^=k; 

numreal=s4 [tempi] [tempj] .r; 
numimag=s4 [tempi ] [tempj] . i; 
fprintf (ofp, ”%e %e\n'’, numreal, numimag) ; 

} /* end for k= . . . */ 

} /* end for j= . . . */ 

} /* end for i= . . . */ 

/* close the output file */ 
fclose (ofp) ; 

} /* end else ... */ 

} /* end if reduce =0 ... */ 

else { /* reduce the amount of output */ 

if ( ( (argc==9) && (argv[7] [1]==' a' ) ) I I ( (argc== 8 ) && (argv[ 6 ] [l]=='a' ) ) ) { 

/* make an ascii output file of the whole spectral plane */ 

/* open the output file and place header information in it */ 
ofp = fopen(outfile,"w"); 
max_num_alf=nprime; 
max_num_f =npr ime ; 

fprintf (ofp, "%i %i\n", max_num_alf ,max_num_f ) ; 
for (i= 0 ; i<max_num_alf ; i++) { 
for (j= 0 ; j<nprime; j++) { 
numreal=s4 [i*halfp] [j] .r; 
numimag=s4 [i*halfp] [j] .i; 

bigmag=sqrt (numreal *numreal + numimag * numimag ) ; 
for (k=l ; k<halfp; k++) { /* look for largest value */ 

tempreal=s4 [i*halfp+k] [j].r; 
tempimag=s4 [i*halfp+k] [j] . i; 

tempmag=sqrt (tempreal*tempreal + tempimag*tempimag) ; 
if (tempmag>bigmag) { /* found a larger value */ 

bigmag=tempmag; 
numreal=tempreal; 
numimag= t emp imag ; 
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} /* end if tempmag>bigmag. . . * 

} /* end if k=l ... */ 

fpnnt f (ofp, " %e %e\n" , numreal, numimag) ; 

} /* end for j=. . . */ 

} /* end for i= . . . */ 

fclose (ofp) ; 

} /* end if ascii... */ 

else { /* make a plot_sxaf reduced output file */ 

/* open the output file and place header information in it * 
ofp = fopen (out file, "w" ) ; 

/* place header information in the file */ 
data_type=2; /* real_or_complex */ 
max_num_alf=npnme; /* num_alpha * / 
alpha_min=0 . 0 ; /* alpha_min */ 

alpha_max=l . 0; /* alpha_max */ 

max_num_f=nprime; /* max_num_f */ 
f_min_all= -0.5; /* f_min_all */ 

f_max_all= 0.5; /* f_max_all */ 

fpnnt f (ofp, " %i %i %e %e %i %e %e\n" , data_type, max_num_aif , alpha_min, 
alpha_max,max_num_f , f_min_all, f_max_all) ; 
for (i=0; i<nprime; i++) { /* group of lines */ 

/* place alpha line header information in file */ 
curr_alf=i; /* alpha */ 
num_f=nprime-i; /* num_f */ 
f_mm= -0 . 5+ (0 . 5*i/npnme) ; /* f_min */ 

f_max= 0 . 5- (0 . 5*i/npnme) ; /* f_max */ 

fprintf (ofp, "%i %i %e %e\n” , curr_alf , num_f , f_min, f_max) ; 
for (k=i; kcnprime; k++) { /* points on the line */ 

tempi= (halfp* (nprime-1) ) - ( (k-i) *halfp) ; 
temp]=k; 

numreal=s4 [tempi ] [tempj] .r; 
numimag=s4 [tempi ] [tempj] . i; 

bigmag=sqrt (numreal '‘'numreal + numimag* numimag) ; 
for (j=l; j<halfp; j++) { /* look for the largest value */ 

tempreal=s4 [tempi+ j ] [tempj] .r; 
tempimag=s4 [tempi+ j ] [tempj] . i; 

tempmag=sqrt (tempreal*tempreal + tempimag*tempimag) ; 
if (tempmag>bigmag) { /* found a larger value */ 

b igmag= t empmag ; 
numreal=tempreal ; 
numimag=t emp imag ; 

} /* end if tempmag> .... */ 

} / * end for j= . . . */ 

fprintf (ofp, " %e %e\n" , numreal, numimag) ; 

} /* end for k=. . . */ 

} /* end for i= . . . */ 

/* close the output file */ 
fclose (ofp) ; 

} /* end else binary... */ 

} /* end reduce the amount of output */ 

} /* end of main */ 
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APPENDIX C. SSCA PROGRAM USE 



Correct commands are: 
for the cross-ssca 

ssca inputfilel inputfile2 outputfile n nprime 1 Oflag 
ssca inputfilel inputfile2 outputfile n nprime 1 Oflag red 

or for the auto- ssca 

ssca inputfilel outputfile n nprime 1 Oflag 
ssca inputfilel outputfile n nprime 1 Oflag red 

where: inputfilel is the file containing the x signal samples 
inputfile2 is the file containing the y signal samples 
outputfile is where the spectrum values will be placed 
n is the number of samples to use from the inputfiles 
and it must be a power of 2 
nprime is the group size of the input datasets and it 
must be a power of 2 

1 is the starting offset of subsequent datasets and 
it must be a power of 2 

Oflag indicates the output filetype, ascii or binary 
-asc for an ascii, MATLAB compatible file, 
full cyclic spectral plane 
-plo for a plot_sxaf , SSPI compatible file, 
half cyclic spectral plane 
red reduces the amount of data output 

note 1: the first two entries in every input file are 

expected to be the datatype and n. datatype = 1 for real 
values, 2 for complex values. n is the number of samples 
contained in the file, one sample per line. 

note 2: error may occur using the reduce option if 

(n-nprime) /nprime is not an integer. 

input file format: 

datatype n 
sample #1 
sample #2 



sample #n 
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output ASCII MATLAB format: 



nprime n 

output (1) (1) 
output (1) (2 ) 



output (nprime) (n) 



output SSPI plot_sxaf format: 

datatype (1 for real, 2 for complex) 
number of alphasalphaminalphamax 
number of f reqsf reqminf reqmax 

value of alpha #1 

number of freqs in #lf reqmin_in_#lf reqmax_in_#l 
output for f reqmin_in_#l 



output for f reqmax_in_#l 



value of alpha #alphamax 

number of freqs in #alphamaxf reqminf reqmax 
output for f reqmin_in_#alphamax 



output for f reqmax_in_#alphamax 
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APPENDIX D. SSCA PROGRAM LISTING 
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^include <stdlib.h> 

^include <staio.h> 

^include <math.h> 

^include "/home3/carter/thesis/SSPI/pam/f ft . c" 

# include " /home3 / carter / thes is / SSP I /pam/ radix . c" 
main (argc, argv) 
int argc; 
char *argv[]; 

{ 

COMPLEX *x, **s3, ♦♦dwnconv; 

COMPLEX *tempfft, *sl, *yl, **s2, ♦window; 
mt l, j, k, type, n, nprime,p, 1, r, direction, norm; 
int a, cross, f ile_n, m, reduce, sizes3, redindex; 
int tempi, temp j, tempk, tempi; 
mt ♦outint; 

float numreal, numimag, convfac, mainfac; 
float bigmag, tempmag, tempreal, tempimag; 
float twopi=6 . 28318530718; 
float *outreal; 
double z,y; 

char *infilel, *infile2, *outfile; 

FILE ♦ifpl, *ifp2, *ofp; 

/* x passes data to/from fft routine 

s3 holds results of operations after first fft 
si holds the x sample values from inputfile 1 

yl holds the y sample values from inputfile 2 

s2 holds the channelized x input data from si 

*/ 

/* 

check for the correct number of input arguments 
the correct commands are : 

ssca inputfile outputfile n nprime 1 Oflag 
ssca inputfile outputfile n nprime 1 Oflag reduce 
ssca inputfilel inputfile2 outputfile n nprime 1 Oflag 

ssca inputfilel mputfile2 outputfile n nprime 1 Oflag reduce 

where: inputfilel is where the x signal samples are expected to be 

inputfile2 is where the y signal samples are expected to be 

outputfile is where the spectrum values will be put 

n is the number of input samples to use 

nprime is the group size of input datasets 
1 is the starting offset of subsequent datasets 
Oflag indicates the output file type 

-asc for an ascii, matlab compatible file 
-plo for an ascii, plot_sxaf compatible file 
reduce is an option to generate a reduced amount of output 

note 1: type and file_n are expected to be on the first line of the mputfile/s 
type = 1 for real, 2 for complex, 

file_n = number of samples following the first line 

note 2: error may occur if (n-npnme) is not evenly divisible by nprime 
when using the reduce option 



maintenance history 

....written by LCDR Nancy J. Carter, USN NPGS Monterey CA 01 October 1992 
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....15 Oct 92.... copied from current version of fam.c to alter to ssca algorithm 
. . .modified to perform ssca per Brown/Gardner/Loomis paper 
....20 Oct 92.... made outputfile ascii format MATLAB compatible 
....22 Oct 92.... made an output format SSPI plot_sxaf compatible 
. . . .24 Nov 92 . . . .added option for a reduced amount of output 
* / 

/* */ 

CHECK INPUT VALUES */ 

/* */ 

if ( (argc 1= 7)&&(argc != 8)&&(argc != 9) ) { 
printf ("ssca . . . fatal error\n" ) ; 

print f(" incorrect number of calling arguments \n" ) ; 

printf (" correct formats are:\n") ; 

printf (" ssca inputfile outputfile n nprime 1 Oflag\n"); 

printf (" ssca inputfile outputfile n nprime 1 Oflag reduce\n"); 

printf (" ssca inputfilel inputfile2 n outputfile nprime 1 Oflag\n"); 

printf (" ssca inputfilel inputfile2 n outputfile nprime 1 Oflag reduce\n") 

printf (" wherein"); 

printf (" inputfilel contains the x signal samplesin"); 

printf (" inputfile2 contains the y signal samplesin"); 

printf (" outputfile will contain the resultsin”); 

printf (" n is the number of input samples to use\n"); 

printf (" nprime is the group size of input datasetsin" ) ; 

printf (" 1 is the offset of subsequent datasetsin" ) ; 

printf (" Oflag is the output file format \n"); 

printf (" Oflag = -asc, a MATLAB file is producedin" ) ; 

printf (" Oflag = -plo, a plot_sxaf file is producedin" ) ; 

printf (" reduce is an option for reduced amount of outputin"); 

exit ( ) ; 

} 

if (argc==7) { /* this is an auto-ssca, no output reduction */ 

cross=0 ; 



inf ilel=argv [1] ; 
outfile=argv[2] ; 
n=atoi (argv [3] ) ; 
nprime=atoi (argv [4] ) ; 
l=atoi (argv [5] ) ; 
reduce=0 ; 

} 

if ( (argc==8) && (argv [7] [0] ! =' r ' ) ) { /* this is a cross-ssca, no data reduction */ 

cross=l; 

inf ilel=argv [ 1 ] ; 
inf ile2=argv [2 ] ; 
outfile=argv [3] ; 
n=atoi (argv [4 ] ) ; 
nprime=atoi (argv [5] ) ; 
l=atoi (argv [ 6] ) ; 
reduce=0; 

if ( (argc==8) && (argv [7 ] [0] ==' r' ) ) { /* this is an auto-ssca, with data reduction */ 

cross=0; 

inf ilel=argv [1 ] ; 
out file=argv [2 ] ; 
n=atoi (argv [3] ) ; 
nprime=atoi (argv (4] ) ; 
l=atoi (argv [ 5 ] ) ; 
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reduce=l; 



if (argc==9) { / * this is a cross-ssca, with data reduction * 

cross=l; 

inf ilel=argv [ 1 ] ; 
inf ile2=argv [2 ] ; 
outf ile=argv [3] ; 
n=atoi (argv [4 ] ) ; 
nprime=atoi (argv [5 ] ) ; 
l=atoi (argv [ 6] ) ; 
reduce=l ; 

} 

/* verify that n, nprime and 1 are powers of 2 



i=n%2 ; 
if (i ! =0 ) { 

printf ("ssca . . . fatal errorXn" ) ; 

printf (" calling argument n is not a power of 2\n"); 



exit ( ) ; 

} 

i=nprime%2 ; 



printf ("ssca. . .fatal errorXn") ; 

printf (" calling argument nprime is not a power of 2\n") ; 



exit ( ) ; 

} 

j-1%2; 
if ( j ! =0) { 



printf ("ssca . . . fatal errorXn" ) ; 

printf (" calling argument 1 is not a power of 2\n"); 

exit () ; 

} 

/* */ 

/* open input file 1, prepare to get the x signal sample data */ 
/* */ 



lfpl = fopen (infilel, "r" ) ; 

fscanf (ifpl, "%i %i", &type, &file_n) ; 

/* verify that file_n is at least equal to n */ 
if ( f ile_n < n) { 

printf ("ssca. . . fatal error\n") ; 

printf (" inputfilel does not contain enough samples\n"); 

fclose (ifpl) ; 
exit ( ) ; 

} 

/* find p - the number of datasets (nprime long) */ 
p= ( (n-npnme) /l) ; 

/* find one dimension of final data array */ 
sizes3=p*l; 

/* */ 

/* READ IN DATA SAMPLES */ 

/* */ 

/* allocate space to read in the x sample values */ 

/* V 

s 1= (COMPLEX*) calloc (n, sizeof (COMPLEX) ) ; 
if ( S 1==NULL) { 

printf ("ssca ... insufficient space to allocate sl\n"); 
exit () ; 



7 ^ 



Dec 9 12:22 1992 strip/ssca.c Page 4 



if (type==l) { * read in the real sample values */ 

for (i=0; i<n; i++) { 

fscanf (ifpl, " %e\n" , Snumreal ) ; 
sl[i] . r=numreal; 
si [i] . 1 = 0 . 0 ; 

} 

} 

else { /* read in the complex sample values */ 

for (i=0; i < n; i++) { 

fscanf (ifpl, " %e %e\n", snumreal, Snumimag) ; 
si [i] .r=numreal; 
sl[i] .i=numimag; 

} 

} 

/* */ 

/* close the input file 41 */ 

*/ 

fclose (ifpl) ; 

/* * / 

if (cross==l) { 

/* */ 

/* open inputfile2, prepare to get the y signal sample data */ 

/ * */ 

ifp2 = fopen (infile2, "r" ) ; 
fscanf (ifp2, " %i %i", &type, &file_n) ; 

/* verify that file_n is at least as big as n */ 
if (file_n < n) { 

printf C’ssca . . . fatal error\n" ) ; 

printfC inputfile2 does not contain enough samples\n") 

fclose (ifp2 ) ; 
exit ( ) ; 

} 

/* */ 

/* allocate space to read in the y sample values */ 

/* */ 

yl= (COMPLEX* ) calloc (n, sizeof (COMPLEX) ) ; 
if (yl==NULL) { 

printf ("ssca ... insufficient space to allocate yl\n" ); 
exit () ; 

} 

/ * */ 

/* read in the complex sample values */ 

V 

for (i=0; i < n; i++) { 

fscanf (ifp2, "%e %e\n" , finumreal, Snumimag) ; 

yl[i] . r=numreal; 
yl[i] .i=numimag; 

} 

/* */ 

/* close the inputfile2 */ 

/* */ 

fclose (ifp2) ; 

} /* end if cross =1 */ 

/* */ 

/* FORM DATASETS */ 
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*/ 

- * allocate space to hold the p-by-nprime datasets 

V 

s2= (COMPLEX**) calloc (p, sizeof (COMPLEX*) ) ; 
if ( s2==NULL) { 

print f ( "ssca ... insufficient space to allocate s2\n"); 
exit ( ) ; 

} 

for (i=0; i<p; i++) { 

s2 [i] = (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 
if (s2 [i] ==NULL) { 

print f ( "ssca ... insufficient space to allocate s2\n"); 
exit ( ) ; 

} 

} 

/* */ 

/* form p datasets of nprime samples from the original sample stream */ 
/* V 

for (i=0; i<p; i++) { 

for ( j=0; j<nprime; j++) { 
k=i*l+ j ; 
s 2 [ i ) [ j ] =s 1 [ Jc ] ; 

} 

} 

/* */ 

/* APPLY WINDOW TO DATASETS */ 

/* */ 

/* allocate space to hold the window multiplicand values */ 

/* */ 

window= (COMPLEX* ) calloc (nprime, sizeof (COMPLEX) ) ; 
if (window==NULL) { 

printf ( "ssca ... insufficient space to allocate sl\n"); 
exit ( ) ; 

} 

/* */ 

/* calculate the window values for the "nprime" sample wide window */ 

/* */ 

for (i=0; Knprime; i++) { 
y= (twopi*i) / (nprime-1) ; 
z=cos (y) ; 

window [i].r= 0.54 - (0.46*z); 
window [i] . i= window [i] . r; 

} 

/* */ 

/* apply the window to rows of s2 containing data */ 

/* */ 

for (i=0; i<p; i++) { 

for (j=0; j<nprime; j++) { 

s2 [i] [ j] .i=s2 [i] [ j] . i*window[ j] . i; 
s2 [ i ] [ j ] . r=s2 [ i ] [ j ] . r *window [ j ] . r ; 

} 

} 

/* */ 

/* FFT EACH DATASET ROW */ 

/* */ 

/* allocate space to hold rows of s2 for passing to the FFT routine */ 
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x= (COMPLEX'"' ) calloc (npnme, sizeof (COMPLEX) ) ; 
if ( x==NULL) { 

printf ("ssca. .. insufficient space to allocate x\n"); 
exit ( ) ; 

*/ 

* allocate space to hold rows of results from the FFT routine */ 
/* */ 

s3= (COMPLEX* * ) calloc (nprime, sizeof (COMPLEX*) ) ; 
if (s3==NULL) { 

print f ("ssca ... insufficient space to allocate s3\n"); 
exit ( ) ; 

) 

for (i=0; i<npnme; i++) { 

s3 [i] = (COMPLEX*) calloc ( (p*l) , sizeof (COMPLEX) ) ; 
if (s 3 [ i ] ==NULL ) { 

printf ("ssca ... insufficient space to allocate s3\n"); 
exit ( ) ; 

} 

} 

/* */ 

/* take an FFT of each row of s2 */ 

/* */ 

/* set values in arguments sent to fft */ 
direction=l ; 
norm=l ; 

for (i=0; i<p; i++) { 

/* copy row of s2 into complex array */ 

for ( j = 0; jcnprime; j++) { 
x [ j]=s2 [i] [ j] ; 

} 

/* go get FFT performed */ 

fft (x, nprime, direction, norm) ; 

/* copy x values into 1th column of s3 */ 

r=i*l; 

for (j=0; jcnprime; j++) { 
s3[j] [r ] =x [ j ] ; 

} 

} 

/* */ 

/* DOWNCONVERT EACH ROW */ 

*/ 

/* allocate space to hold the downconversion multiplicands */ 

/* */ 

awnconv= (COMPLEX**) calloc (p, sizeof (COMPLEX*) ) ; 
if (dwnconv==NULL) { 

print f (" ssca ... insufficient space to allocate dwnconv\n"); 
exit () ; 

} 

for (i=0; i<p; i++) { 

dwnconv [i] = (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ) ; 
if (dwnconv [i] ==NULL) { 

print f ("ssca ... insufficient space to allocate dwnconv\n"); 
exit ( ) ; 

} 
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/* downconvert each of the transform sequences *■ 

. * calculate the down conversion multipliers * 

/ * * / 
mainfac=twopi*l/nprime; 
for (m=0; m<p; m++) { 
for (k=0; k<nprime; k++) { 

convfac=mainfac* (k- (nprime/2) +1) *m; 
dwnconv[m] [k] . r=cos (convfac) ; 
dwnconv[m] [k] . i= (-1 . 0) *sin (convfac) ; 

} 

} 

/* multiply downconversion factor against each frequency value */ 
for (1=0; i<p; i++) { 
tempi=i*l; 

for ( j=0; j<nprime; ]++) { 

numreal=s3 [ j ] [tempi] . r *dwnconv [ i ] [j] .r - s 3 [ j ] [tempi] . i*dwnconv [ i] [j] . 1 ; 

s3 [ j ] [tempi ] . i=s3 [ j ] [tempi ] . r ‘dwnconv [ i] [ j ] . i + s3 [ j ] [tempi] .i T dwnconv[i] [ j J . r 

s3[j] [tempi] .r=numreal; 

} 

} 

/* */ 

/* REPLICATE COLUMNS */ 

/* */ 

for (i=0; i<p; i++) { 
r=i*l; 

/* copy column into 1-1 adjacent columns */ 
for (j=0; jcnprime; j++) { 

for (k=l ; k<l; k++) { 
s3[ j] [r+k] =s3 [ j ] [r] ; 

} /* end for k= . . . */ 

} / * end for j= . . . * / 

} /* end for i= . . . */ 

/* */ 

/* MULTIPLY COLUMNS */ 

/* */ 

if (cross==0) { /* auto-ssca */ 

/* */ 

/* multiply columns of s3 with conjugate value of si */ 

/* */ 

for (a=0; a<sizes3; a++) { 
for (j=0; j<nprime; j++) { 

s 3 t j ] [a] . r=(s3 [ j] [a] .r*sl [a] .r) + (s3 [ j] [a] . i*sl [a] . i) ; 
s 3 [ j ] [a] . i= (s3 [ j ] [a] . i*sl [a] . r)- (s3 [ j] [a] . r*sl [a] . i) ; 

} 

} /* end for a=0 ... */ 

} /* end auto-ssca */ 

else { /* cross-ssca */ 

/* */ 

/* multiply columns of s3 with conjugate value of yl */ 

/* */ 

for (a=0; a<sizes3; a++) { 
for ( j=0; j<nprime; j++) { 

s3 [ j ] [a] . r= (s3 [ j ] [a] .r*yl[a] .r)+(s3[j] [a] . i*yl[a] .i); 
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s3 [ j] [a] . i= (s3 [ j ] [a] . i*yl [a] . r)-(s3 [ j] [a] .r*yl [a] . i) ; 

} 

} / * end for a=0 ... * / 

/* end cross-ssca */ 

* */ 

/* FFT EACH ROW */ 

/* */ 

* */ 

- * allocate space to hold rows of s3 for passing to the FFT routine */ 

*/ 

tempf ft= (COMPLEX* ) calloc (sizes3, sizeof (COMPLEX) ) ; 
if (tempf ft ==NULL) { 

printf ("ssca. .. insufficient space to allocate tempf ft %i by l\n" / sizes3 
exit ( ) ; 

} 

/* set values in arguments sent to fft */ 

direction=l ; 

norm=l; 

for (j = 0; j<npnme; j++) { 

/* copy row of s3 into tempf ft */ 
for (k=0; k<sizes3; k++) { 
tempf ft [k] =s3 [ j ] [k] ; 

} 

/* go get FFT performed */ 

fft (tempf ft, sizes3, direction, norm) ; 

/* copy tempfft result values back into a row of s3 */ 

for (k=0; k<sizes3; k++) { 
s3[j] [k] =tempf ft [k] ; 

} 



} /* end for j=0 . . . 


*/ 


/* 


*/ 


/ * OUTPUT 


*/ 


/* 


*/ 


if (reduce==0) { /* 


no output reduction 



if ( ( (argc==7)&& (argv[6] [l]=='a' ) ) I I ( (argc==6) && (argv [5] [ 1 ] ==' a' ) ) ) { 

/* make an ascii output file full spectral plane */ 

/* open the output file and place header information in it */ 

ofp = fopen (out file, "w" ) ; 

fprintf (ofp, "%i %i\n", nprime, sizes3) ; 

/* write out data values to file for plotting */ 
for (i=0; Knprime; i++) { 
for ( j=0; j<sizes3; j++) { 

fprintf (ofp, " %e %e\n", s3 [i] [ j ] . r, s3 [i] [ j ] . i) ; 

} /* for j=0 . . . */ 

} /* for i=0 ... * / 

/* close the output file */ 
fclose (ofp) ; 

} /* if ascii ... */ 

else { /* make a plot_sxaf output file half spectral plane */ 

/* open the output file and place header information in it */ 
ofp = fopen (outfile, "w" ) ; 
a=0; /* alpha = 0 */ 

fprintf (ofp, "%i %i %e %e %i %e %e\n", 2, sizes3, 0 . 0, 1 . 0, nprime, -0 . 5, 0 . 5) 
for (i=0; i<nprime; i++) { /* group of lines */ 

for (j=0; j< (sizes3/nprime) / j++) { /* line in the group */ 

tempreal= -0 . 5+ (0 . 5*i/nprime) ; 
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tempimag= 0 . 5- ( 0 . 5*i/npnme) ; 

fprint f (ofp, " %i %i %e %e\n" , a, (nprime-i) , tempreal, tempimag) ; 
a=a+l ; 

for (k= 0 ; k< (npnme-i) ; k++) { /* points on the line */ 

tempi=nprime-k-l ; 
temp 3 = (i+k) *sizes3/npnme + :j ; 
numreal=s3 [tempi] [tempj] .r; 
numimag=s3 [tempi] [tempj] . i; 
fpnntf (ofp, "%e %e\n" , numreal, numimag) ; 

} /* end for k=. . . */ 

} /* end for 3 =. . . */ 

} /* end for i=. . . */ 

/* close the output file */ 
fclose (ofp) ; 

} /* end if_else */ 

} /* end if reduce= 0 . . . */ 

if (reduce==l) { /* reduce the amount of output data */ 

if ( ( (argc== 8 ) && (argv [ 6 ] [1] ==' a' ) ) I I ( (argc==7) && (argv [5] [1] ==' a' ) ) ) { 

/* make an ascii output file full spectral plane */ 

/* open the output file and place header information m it */ 
ofp = f open (out file, "w" ) ; 
fpnntf (ofp, "%i %i\n" , npnme, p) ; 
redindex=l; 

/* write out data values to file for plotting */ 
for (i= 0 ; i<nprime; i++) { 
for ( j= 0 ; j<p; j++) { 

numreal=s3 [i] [j*redindex] .r; 
numimag=s3 [i] [j*redindex] . i; 

bigmag=sqrt (numreal*numreal + numimag* numimag) ; 
for (k=l; kcredindex; k++) { /* look for largest value */ 

tempreal=s3 [i] [j*redindex +k] .r; 
tempimag=s3 [i] [j*redindex +k] . i; 

tempmag=sqrt (tempreal* tempreal + tempimag* tempimag) ; 
if (tempmag>bigmag) { /* found a bigger value, swap */ 

b igmag=t empmag ; 
numreal=tempreal ; 
numimag=t emp imag ; 

} /* end if tempmag>.... */ 

} /* end for k=l... */ 

fpnntf (ofp, "%e %e\n", numreal, numimag) ; 

} /* for j= 0 . . . */ 

} /* for i= 0 ... */ 

/* close the output file */ 
fclose (ofp) ; 

} /* if ascii ... */ 

else { /* make a plot_sxaf output file with reduction half spectral plane */ 

/* open the output file and place header information in it */ 
ofp = fopen (outf ile, "w" ) ; 
redindex=sizes3/npnme; 
a= 0 ; /* alpha = 0 */ 

fprintf (ofp, "%i %i %e %e %i %e %e\n", 2 , nprime, 0 . 0 , 1 . 0 , npnme, - 0 . 5, 0 . 5) ; 
for (i= 0 ; icnprime; i++) { /* nprime alpha levels */ 

tempreal= - 0 . 5+ ( 0 . 5*i/npnme) ; 
tempimag= 0 . 5- ( 0 . 5 *i/nprime) ; 

fprintf (ofp, "%i %i %e %e\n", a, (nprime-i) , tempreal, tempimag) ; 
a=a+l; 
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for (k=0; k< (nprime-i ) ; k++) { - ' points at the alpha level 

t emp i =np r ime-k- 1 ; 
temp j= (i+k) *redindex; 
numreal=s3 [tempi] [temp]] . r; 
numimag=s3 [tempi ] [temp]] . i; 

bigmag=sqrt (numreal*numreal + numimag*numimag) ; 
for (j=l; j<redindex; j++) { /* line in the group */ 

tempreal=s3 [tempi ] [temp] +j].r; 
tempimag=s3 [tempi ] [tempj +j] . i; 

tempmag=sqrt (tempreal*tempreal + tempimag^tempimag) ; 
if (tempmag>bigmag) { /* found a bigger value, swap */ 

b igmag=t empmag ; 
numreal=tempreal ; 
numimag=t empimag; 

} /* end if tempmag>.... */ 

} /* end for j=. . . */ 

fprint f (ofp, "%e %e\n" , numreal, numimag) ; 

} /* end for k=. . . + / 

} /* end for i=. . . */ 

/* close the output file */ 
fclose (ofp) ; 

} /* end if_else */ 

} /* end if reduce=l... */ 

} 
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APPENDIX E. SUBFAM PROGRAM USE 



Correct command is : 

subfam n inputfilel offsetl freqzerol conjl inputfile2 offset2 
freqzero2 conj2 outputfile IOflag Itypes 

where : 

n is the number of samples to use from each input file, 
it must be a power of 2 

inputfilel is the first file containing signal samples 

offsetl is the initial sample position offset location 
in inputfilel 

freqzerol is the flag indicating which frequency 
modification to perform on inputfilel, options are: 
zn - zero neg freqs, shift pos freqs down 
zp - zero pos freqs, shift neg freqs up 
zz - don't do anything to frequency components 

conjl is the flag indicating if conjugation of the 
data from inputfilel is to be performed prior to 
multiplication with data from inputfile2, options: 
y - yes, conjugate data 
n - do not conjugate data 

inputfile2 is the second file containing signal 
samples 

offset2 is the initial sample position offset location 
in inputfile2 

freqzero2 is the flag indicating which frequency 
modification to perform on inputfile2, options are: 
zn - zero neg freqs, shift pos freqs down 
zp - zero pos freqs, shift neg freqs up 
zz - don't do anything to frequency components 

conj2 is the flag indicating if conjugation of the 
data from inputfile2 is to be performed prior to 
multiplication with data from inputfilel, options: 
y - yes, conjugate data 
n - do not conjugate data 

outputfile is where the spectrum values will be placed 
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IOflag indicates the datatypes of inputfilel, 

inputfile2 and the outputfile. three characters are 
given to indicate the datatype of each file, 
valid datatypes are: 
a - file is of type ascii 
b - file is of type binary 

ex. "abb" means inputfilel is of type ascii, 
inputfile2 is of type binary and the outputfile 
is to be of type binary 

Itypes indicates the types of data inside inputfilel 
and inputfile2. the outputfile is always complex 
floating point, real and imaginary components on the 
same line. two letters must be given to indicate 
the datatype for each input file. all types are 
converted internally to complex floating point, 
valid types are: 

r - single real floating point samples per line 
i - single integer sample per line 
c - two real floating point numbers per line, 
representing the real and imaginary parts 
of a single sample 

ex. "ic" means inputfilel has integer samples 
and inputfile2 has complex floating point samples 



subfam 8192 datal.dat 64 zn n data2.dat 2156 zp y 
outputA.dat bbb rr 



example : 



input file format: 



output file format: 

magnitude #1 phase #1 
magnitude #2 phase #2 



sample #1 
sample #2 



sample #n 



magnitude #n phase #n 
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^include <stdlib.h> 

^include <stdio.h> 

# include <math.h> 

$ include " /home3/carter/thesis/SSPI/pam/f ft . c" 

# include " /home3 /carter /thesis /SSP I /pam/radix. c" 

/* values used in up/down conversion */ 

^define INTFREQ 49375000.0 /* intermediate signal freq in Hz */ 

^define SAMPTIME 0.0000051282 /* ong signal sampling period in sec */ 

main (argc, argv) 
int argc; 
char *argv [ ] ; 

{ 

COMPLEX * s 1 , * s2 , *s3; 

int i, n, direction, norm, of f 1, of f2, quarter_n, scaninteger, *readinteger, *z; 
float tempreal, tempimag, outputmag, outputphase, numreal, intfreq, t, convfac; 
float *readreal, *readimag, *y; 

char *infilel, *infile2, *outfile, *conjl, *conj2, *freqzerol, *freqzero2, *ioflag, *itypes; 
FILE *ifpl, *ifp2, *ofp; 

/* 

subfam is a program which uses signal samples from two input files to calculate 
spectral values in a single diamond of the frequency-alpha plane 

the command is : 

subfam n inputfilel offsetl freqzerol conjl inputfile2 offset2 freqzero2 
conj2 outputfile IOflag Itypes 

n is the number of samples to use from each inputfile, must be pwr of 2 

inputfilel is the name of the first input file containing signal samples 
offsetl is the initial sample position offset location in inputfilel 
freqzerol is the flag indicating which frequency modification to perform 
on inputfilel, valid options are: 

zn - zero negative frequency components, shift positive 
components down 

zp - zero positive frequency components, shift negative 
components up 

zz - don't do anything to frequency components 
conjl is the flag indicating if conjugation of the data from inputfilel is 
to be performed prior to multiplication with data from mputfile2. 
valid options are: 
y - yes conjugate data 
n - do not conjugate data 

mputfile2 is the name of the second input file containing signal samples 
offset2 is the initial sample position offset location in inputfile2 
freqzero2 is the flag indicating which frequency modification to perform 
on input file2, valid options are: 

zn - zero negative frequency components, shift positive 
components down 

zp - zero positive frequency components, shift negative 
components up 

zz - don' t do anything to frequency components 
conj2 is the flag indicating if conjugation of the data from inputfile2 is 
to be performed prior to multiplication with data from inputfilel. 
valid options are: 
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y - yes conjugate data 
n - do not conjugate data 

outputfile is the filename where the spectrum values will be placed 

IOflag indicates the datatypes of inputfilel, inputfile2 and outputfile. 
three letters must be given to indicate the datatypes for all files, 
valid datatypes are: 
a - file is of type ascii 
b - file is of type binary 

ex. "abb" means inputfilel is of type ascii, 
inputfile2 is of type binary, and 
the outputfile is to be of type binary 

Itypes indicates the number types of data inside inputfilel and 

inputfile2. the outputfile is always complex floating point, 
real and imaginary components on the same line. two letters 
must be given to indicate the number type for each input file, 
all number types are converted to complex floating point, 
valid number types are: 

r - single real floating point sample per line 
i - single integer sample per line 

c - two floating point numbers per sample on each line 
ex. "ic" means inputfilel has integer samples and 

inputfile2 has complex floating point samples 



sample useage: 

subfam 8192 datal.dat 64 zn n data2.dat 2156 zp y outputA.dat bbb rr 
maintenance history 

. . . created 01 October 1992, LCDR Nancy J. Carter, USN at NPGS Monterey CA 

... 10.15.92 - added freqzero options - njc 

... 10.16.92 - added error check on conjl and conj2 - njc 

... 10.25.92 - changed down/up shifting to same code as in old version 



V 

/* check for the correct number of input arguments */ 

if (argc 1= 13) { 

printf ( "subfam. .. fatal error ... incorrect number of calling argument s \n ") ; 
print f ( "\n" ) ; 

printf (".... correct format is: \n"); 

printf (" subfam n inputfilel offsetl freqzerol conjl inputfile2 offset2 \n"); 
printf (" freqzero2 conj2 outputfile IOflag Itypes\n"); 

printf ( "\n" ) ; 
printf ("\n") ; 

printf (" n is the number of samples to use from each inputfile, must be pwr of 2\n" 

printf (" inputfilel is a file containing the signal samples\n"); 

printf (" offsetl is the initial sample position offset in input filel\n" ) ; 

printf (" freqzerol indicates type of freq zeroing and shifting to do on inputfilelX 

printf (" zn....zero negative freqs and shift down\n"); 

printf (" zp....zero positive freqs and shift up\n"); 

printf (" zz....do not zero anything, do not shift\n"); 

printf (" conjl indicates if conjugation of data from inputfilel is\n"); 

printf (" desired before multiplication^"); 

printf (" y....yes conjugate data\n"); 
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printfC n no do not conjugate dataxn"); 

printfC usual setup is conjl = n, conj2 = y\n"); 

printfC inputfile2 is a file containing the signal samplesXn"); 

printfC' offset2 is the initial sample position offset in input file2\n" ) ; 

printfC freqzero2 indicates type of freq zeroing and shifting to do on inputfile 

printfC* zn zero negative freqs and shift downXn"); 

printfC* zp....zero positive freqs and shift up\n" ); 

printfC’ zz....do not zero anything, do not shiftXn"); 

printfC' conj2 indicates if conjugation of data from mputfile2 is\n"); 

printfC desired before multiplicationXn" ) ; 

printfC' y....yes conjugate dataXn"); 

printfC n no do not conjugate dataXn"); 

printfC* usual setup is conjl = n, conj2 = y\n"); 

printfC' outputfile is where the spectrum values will be placedXn"); 

printfC* IOflag indicates datatypes of input filel, input file2 and output file\n" ) ; 

pnntf (" aaa ascii, ascii, asciiXn") ; 

printf ( " aab ascii, ascii, binaryXn" ) ; 

pnntf ( " aba ascii, binary, asciiXn" ) ; 

printf (" abb ascii, binary, binaryXn" ) ; 

printf ( " baa ascii, ascii, as cii\n" ) ; 

printf (" bab ascii, ascii, binaryXn" ) ; 

printf (" bba binary, binary, asciiXn" ) ; 

printf (" bbb binary, binary, binary \n" ) ; 

printfC' Itypes indicates number types of input samplesXn"); 

printfC* r single real float\n"); 

printfC' i single integerXn"); 

printfC' c complex, two floatsXn"); 

printf ("\n" ) ; 

print f ( " . . . . sample useage : \n" ) ; 

printfC' subfam 8192 datal.dat 64 zn n data2.dat 2156 zp y outputA.dat bbb rr\n") ; 
exit ( ) ; 

} 

n=atoi (argv [1] ) ; 
inf ilel=argv [2] ; 
off l=atoi (argv [3] ) ; 
freqzerol=argv [4] ; 
conjl=argv [5] ; 
inf ile2=argv [ 6] ; 
of f2=atoi (argv [7 ] ) ; 
freqzero2=argv[8] ; 
conj2=argv [9] ; 
outfile=argv [10] ; 
io f lag=argv [11] ; 
itypes=argv [ 12 ] ; 






/* verify that n is a power of 2 
i=n%2 ; 
if (i ! =0) { 

printf ( "subfam. . . fatal error\n" ) ; 

printfC' calling argument n is not a power of 2\n") ; 

exit ( ) ; 



} 



verify that conjl & conj2 are y or n */ 
if ( ( (conjl [0] ! =' y' ) && (conjl [0] 1=' n' ) ) I | ( (conj2 [0] \=' y' ) && (conj2 [0] !=' n' ) ) ) { 

printf ("subfam. . . fatal error\n") ; 

printfC conjl and conj2 must be set to y or n\n"); 
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exit () ; 

} 

/* verify that itypes are r, i or c * > 

if ( ( (itypes [0] !=' r' ) && (itypes [0] ! =' i' ) && (itypes [0] c' )) I 
( (itypes [1] !=' r' ) && (itypes [1] !=' i' ) && (itypes [1] !=' c' ) ) ) { 
pnntf ("subfam. . . fatal error\n") ; 



printf (" Itypes must be set to r, i or c\n") ; 

exit () ; 

} 

/* */ 

/* allocate space for the 1-D data arrays */ 

/* V 



s 1= (COMPLEX* ) calloc (n, sizeof (COMPLEX) ) ; 
if (sl==NULL) { 

pnntf ("subfam. .. insufficient space to allocate sl\n"); 
exit () ; 

} 

s2= (COMPLEX* ) calloc (n, sizeof (COMPLEX) ) ; 
if (s2==NULL) { 

printf ("subfam. .. insufficient space to allocate s2\n" ); 
exit ( ) ; 

} 

/* allocate space to read in binary data */ 

y= (float *)malloc (1* sizeof (float) ) ; 
z= (int* )malloc (l*sizeof (int ) ) ; 

/* */ 

/* INPUT */ 

/* */ 

/* open input files, get the signal samples */ 

/* */ 

switch (ioflag[0]) { 

case ' a' : /* inputfilel is an ASCII file */ 

if (( ifpl = fopen (infilel, "r" ) ) == NULL ) { 

printf ("subfam. . .unable to open ascii inputfilel\n" ) / 
exit () ; 

} 

/* move to the desired starting offset position in the file */ 
for (i=0; i<of fl; i++) { 
switch (itypes [0]) { 

case ' r' : /* single real float per line */ 

if (fscanf (ifpl, "%e", Stempreal) 1=1) { 

printf ("subfam. . .EOF encountered while attempting to reach offsetl\n"); 
fclose (ifpl) ; 
exit () ; 

} 

break; 

case ' i' : /* single integer per line */ 

if (fscanf (ifpl, "%i", scaninteger) 1=1) { 

print f ( "subfam. . .EOF encountered while attempting to reach offsetl\n"); 
fclose (ifpl) ; 
exit () ; 

} 

break; 

default: /* complex, two floats per line */ 

if (fscanf (ifpl', "%e %e" , stempreal, &tempimag) 1 =2) { 

printf ("subfam. . .EOF encountered while attempting to reach offsetl\n"); 
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fclose (ifpl) ; 
exit () ; 

} 

} /* end switch itypes[0] */ 

} /* end if i=0;i<offl */ 

* read in n data samples */ 
for (i=0; i<n; i++) { 
switch (itypes[0]) { 

case ' r' : /* single real float per line */ 

if (fscanf (ifpl, " %e", &sl [i] . r) 1=1) { 

pnntf ("subfam. . .EOF reached while reading mputf ilel\n" ) ; 
fclose (ifpl) ; 
exit ( ) ; 

} 

break; 

case ' i' : /* single integer per line */ 

if ( fscanf (ifpl, "%i" , scaninteger ) I =1) { 

printf ("subfam. . .EOF reached while reading input filel\n" ) ; 
fclose (ifpl) ; 
exit () ; 

} 

si [i] . r=scaninteger; 
break; 

default: /* complex, two floats per line */ 

if ( fscanf (ifpl, " %e %e" , &sl [ i] . r, &sl [i] . i) ! =2) { 

printf ("subfam. . .EOF reached while reading inputf ilelXn" ) ; 
fclose (ifpl) ; 
exit () ; 

} 

} /* end switch itypes[0] */ 

} /* end for i=0 i<n */ 

fclose (ifpl ) ; 
break; 

case 'b' : /* inputf ilel is a BINARY file */ 

if (( ifpl = fopen (infilel, "rb" ) ) == NULL ) { 

printf ("subfam. . .unable to open binary input filel\n" ) ; 
exit () ; 

} 

/* move to the desired starting offset position in the file */ 
for (i=0; i<of f 1; i++) { 
switch (itypesfO]) { 

case ' r' : /* single real float per line */ 

if (fread (y,sizeof(*y),l,ifpl) 1=1) { 

printf ("subfam. . .EOF encountered while attempting to reach offsetl\n") 
fclose (ifpl) ; 
exit () ; 

} 

break; 

case ' i' : /* single integer per line */ 

if (fread (z, sizeof (*z) , 1, ifpl) 1=1) { 

printf ("subfam. . .EOF encountered while attempting to reach offsetlXn") 
fclose (ifpl) ; 
exit () ; 

} 
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break; 

default: / * complex, two floats per line T 

if (fread(y, sizeof (*y) , 1, ifpl) 1=1) { 

pnntf ( "subfam. . .EOF encountered while attempting to reach offsetl\n") 
fclose (ifpl) ; 
exit ( ) ; 

} 

if (fread(y, sizeof (*y) , 1, ifpl) 1=1) { 

pnntf ("subfam. . .EOF encountered while attempting to reach offsetl\n") 
fclose (ifpl) ; 
exit ( ) ; 

} 

} /* end switch itypes[0] * / 

} /* end if i=0;i<offl */ 

/* read m n data samples */ 
for (i=0; i<n; i++) { 
switch (itypes[0]) { 

case ' r' : /* single real float per line */ 

if ( fread (y, sizeof ( *y) , 1, ifpl) ! =1) { 

printf ("subfam. . .EOF reached while reading input filel\n" ) ; 
fclose (ifpl) ; 
exit ( ) ; 

} 

si [i] . r=*y ; 
break; 

case ' i' : /* single integer per line */ 

if ( fread (z, sizeof (*z) , 1, ifpl) 1=1) { 

printf ("subfam. . .EOF reached while reading inputfilelVn" ) ; 
fclose (ifpl) ; 
exit ( ) ; 

} 

si [i] . r=*z; 
break; 

default: /* complex, two floats per line */ 

if ( fread (y, sizeof(*y),l,ifpl) 1=1) { 

printf ("subfam. . .EOF reached while reading input filel\n" ) ; 
fclose (ifpl) ; 
exit () ; 

} 

si [i] . r=*y; 

if ( fread (y, sizeof (*y) , 1, ifpl) 1=1) { 

printf ("subfam. . .EOF reached while reading inputf ilel\n" ) ; 
fclose (ifpl) ; 
exit () ; 

} 

si [i] . r=*y ; 

} /* end switch itypes[0] */ 

} /* end for i=0 i<n */ 

fclose (ifpl) ; 
break; 

default: /* invalid format specified for inputfilel */ 

printf ("subfam. .. invalid format specified for inputf ilel\n" ) ; 
exit ( ) ; 

} /* end switch ioflag[0] */ 
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switch (ioflag[l]) { 

case 'a ' : /* inputfile2 is an ASCII file */ 

if (( ifp2 = fopen (inf ile2, "r" ) ) == NULL ) { 

printf ("subfam. . .unable to open ascii input file2\n" ) ; 
exit ( ) ; 

} 

/* move to the desired starting offset position in the file *' 
for (i=0; i<of f 2 ; i++) { 
switch (itypes[l]) { 

case ' r' : /* single real float per line */ 

if ( fscanf (ifp2, " %e" , Stempreal) t =1) { 

printf ("subfam. . .EOF encountered while attempting to reach offset2\n") 
fclose (ifp2 ) ; 
exit () ; 

} 

break; 

case ' i' : /* single integer per line */ 

if (fscanf (ifp2, "%i", scaninteger) !=1) { 

printf ("subfam. . .EOF encountered while attempting to reach offset2\n") 
fclose (ifp2 ) ; 
exit () ; 

} 

break; 

default: /* complex, two floats per line */ 

if (fscanf (ifp2, "%e %e", Stempreal, Stempimag) !=2) { 

printf ("subfam. . .EOF encountered while attempting to reach offset2\n") 
fclose (ifp2 ) ; 
exit () ; 

} 

} /* end switch itypes[l] */ 

} /* end if i=0;i<off2 */ 

/* read in n data samples */ 
for (i=0 ; i<n; i++) { 
switch (itypes[l]) { 

case ' r' : /* single real float per line */ 

if (fscanf (ifp2, "%e", &s2 [i] .r) 1=1) { 

printf ("subfam. . .EOF reached while reading inputfile2\n" ) ; 
fclose (ifp2) ; 
exit ( ) ; 

} 

break; 

case ' i' : /* single integer per line */ 

if (fscanf (ifp2, " %i" , scaninteger) ! =1) { 

printf ( "subfam. . .EOF reached while reading inputf ile2\n" ) ; 
fclose (ifp2) ; 
exit () ; 

} 

s 2 £ i ] . r=scaninteger ; 

break; 

default: /* complex, two floats per line */ 

if (fscanf (ifp2, "%e %e", &s2 [i] . r, &s2 [i] .i) 1=2) { 

printf ("subfam. . .EOF reached while reading inputf ile2\n" ) ; 
fclose (ifp2 ) ; 
exit ( ) ; 

} 
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} /* end switch itypes[l] V 

} /* end for i=0 Kn */ 

f close (ifp2) ; 
break; 

case 'b' : /* inputfile2 is a BINARY file */ 

if (( ifp2 = fopen (infile2, "rb") ) == NULL ) { 

print f ("subf am. . .unable to open binary input file2\n" ) ; 
exit ( ) ; 

} 

/* move to the desired starting offset position in the file */ 
for (i=0; i<of f2; i++) { 
switch (itypes[l]) { 

case ' r' : /* single real float per line */ 

if (fread (readreal, sizeof ( *readreal) , 1, ifp2) 1=1) { 

print f ( "subfam. . .EOF encountered while attempting to reach offset2\n"); 
fclose (ifp2) ; 
exit ( ) ; 

} 

break; 

case ' i' : /* single integer per line */ 

if ( fread (readinteger , sizeof ( *readinteger) , 1, ifp2) !=1) { 

printf ( "subfam. . .EOF encountered while attempting to reach offset2\n"); 
fclose (ifp2) ; 
exit ( ) ; 

} 

break; 

default: /* complex, two floats per line */ 

if (fread (readreal, sizeof (*readreal) , 1, ifp2) !=1) { 

printf ("subfam. . .EOF encountered while attempting to reach offset2\n"); 
fclose (ifp2) ; 
exit ( ) ; 

} 

if (fread (readimag, sizeof (*readimag) , 1, ifp2) !=1) { 

printf (" subfam. . .EOF encountered while attempting to reach of fset2\n" ) ; 
fclose (ifp2) ; 
exit ( ) ; 

} 

} /* end switch itypes[l] */ 

} /* end if i=0;i<off2 */ 

/* read in n data samples */ 
for (i=0; i<n; i++) { 
switch (itypes[l]) { 

case ' r' : /* single real float per line */ 

if ( fread (readreal, sizeof ( * readreal) , 1, ifp2) ! =1) { 

printf ( "subfam. . .EOF reached while reading inputf ile2\n" ) ; 
fclose (ifp2) ; 
exit ( ) ; 

} 

s2 [i] .r=*readreal; 
break; 

case ' i' : /* single integer per line */ 

if ( fread (readinteger, sizeof ( *readinteger) , 1, ifp2) 1=1) { 

printf ("subfam. . .EOF reached while reading input file2\n" ) ; 
fclose (ifp2) ; 
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exit 0 ; 

} 

scaninteger=*readinteger ; 
s2[i] . r=scanmteger ; 
break; 

default: /* complex, two floats per line */ 

if (fread (readreal, sizeof (*readreal) , 1, ifp2) !=1) { 

print f ("subf am. . .EOF reached while reading inputf ile2\n" ) 
fclose (ifp2 ) ; 
exit 0 ; 

} 

s2[i] . r=*readreal; 

if (fread (readimag, sizeof (*readimag) , 1, ifp2) ! =1) { 

print f ( "subf am. . .EOF reached while reading inputf ile2\n" ) 
fclose (ifp2) ; 
exit () ; 

} 

s2[i] . i=*readimag; 

} /* end switch itypes[l] */ 

} /* end for i=0 i<n */ 

fclose (ifp2 ) ; 
break; 

default: /* invalid format specified for inputfile2 */ 

printf ("subf am. .. invalid format specified for input file2\n" ) ; 
exit () ; 

, /* end switch ioflag[l] */ 

/ * * / 

/* ZEROIZE REQUESTED FREQUENCY SIDE AND SHIFT APPROPRIATELY */ 

*/ 



quarter_n=n/4; 

switch (freqzerol [1] ) { 

case ' n' : /* zero negative and shift down */ 

direction=l; /* request a forward transform time-to-freq */ 
norm=0; 

f ft (si, n, direction, norm) ; /* results returned in si */ 

for (i=0; i<(n/2); i++) { /* zero lower half */ 
si [i] . r=0; 
si [i] . i=0; 

} 

direction=-l; /* request an inverse transform freq-to-time */ 
norm=2 ; 

f ft (si, n, direction, norm) ; /* results returned in si */ 

* downconvert */ 

t=0.0; /* starting time */ 

for (i=0; i<n; i++) { 

convfac=TWOPI* INTFREQ*t ; 

tempreal=cos (convfac) ; 

tempimag= (-1) *sin (convfac) ; 

numreal=sl [i] . r*tempreal - si [i] . i*tempimag; 

si [i] . i=sl [i] . r*tempimag+sl [i] . i*tempreal; 

s 1 [ i ] .r=numreal; 

t =t+ SAMP TIME; 

} 

break; 

case 'p' : /* zero positive and shift up */ 

direction=l; /* request a forward transform time-to-freq */ 
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norm=0; 

f ft (si, n, direction, norm) ; /* results returned m si V 

for (i= (n/2 ) ; i<n; i++) { /* zero upper half */ 
si [i] . r=0; 
si [i] . i=0; 

} 

direct ion=-l ; /* request an inverse transform freq-to-time */ 

norm=2 ; 

f ft (si, n, direction, norm) ; /* results returned in si V 
/* upconvert */ 

t=0.0; /* starting time */ 

for (i=0; i<n; i++) { 

convf ac=TWOP I * INTFREQ*t ; 
tempreal=cos (convfac) ; 
tempimag=sin (convfac) ; 

numreal=sl [i] . r+tempreal - si [i] . l+tempimag; 
si [i] . i=sl [i] . r*tempimag+sl [i] . i*tempreal; 
s 1 [ i ] . r=numreal; 
t=t+S AMPTIME ; 

} 

break; 

case ' z' : /* do nothing */ 

break; 

default: /* invalid format specified for freqzerol */ 

printf ("subfam. .. invalid format specified for freqzerol\n" ) ; 
exit () ; 

} /* end switch */ 
switch (freqzero2 [1] ) { 

case 'n' : /* zero negative and shift down */ 

direction=l; /* request a forward transform time-to-freq */ 
norm=0 ; 

f ft (s2, n, direction, norm) ; /* results returned in s2 */ 

for (i=0; i<(n/2); i++) { /* zero lower half */ 

s2 [i] . r=0; 
s2 [i] . i=0; 

} 

direct ion=-l; /* request an inverse transform freq-to-time */ 
norm=2; 

fft (s2, n, direction, norm) ; /* results returned in s2 */ 

/* downconvert */ 

t=0.0; /* starting time */ 

for (i=0; i<n; i++) { 

convf ac=TWOP I * INTFREQ *t ; 

tempreal=cos (convfac) ; 

tempimag== (-1) *sin (convfac) ; 

numreal=s2 [i] . r*tempreal - s2 [i] . i*tempimag; 

s2 [i] .i=s2 [i] . r*tempimag+s2 [i] . i*tempreal; 

s2 [i] .r=numreal; 

t=t+SAMPTIME; 

} 

break; 

case ' p' : /* zero positive and shift up */ 

direction=l; /* request a forward transform time-to-freq */ 
norm=0 ; 

fft (s2 , n, direct ion, norm) ; /* results returned in s2 */ 

for (i=(n/2); i<n; i++) { /* zero upper half */ 
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s2 [i] . r=0 ; 
s2 [i] .i=0; 

} 

direct. ion=- 1 ; /* request an inverse transform freq-to-time 

norm=2 ; 

f ft (s2, n, direction, norm) ; /* results returned in s2 */ 
upconvert */ 

t=0.0; /* starting time */ 

for (i=0; i<n; i++) { 

convf ac=TWOP I * INTFREQ * t ; 
tempreal=cos (convfac) ; 
tempimag=sin (convfac) ; 

numreal=s2 [i] . r*tempreal - s2 [ i ] . i*tempimag; 
s2 [i] . i=s2 [i] . r *tempimag+s2 [i] . i*tempreal; 
s2 [i] . r=numreal; 
t=t+SAMPTIME; 

} 

break; 

case ' z' : /* do nothing */ 

break; 

default: /* invalid format specified for freqzero2 */ 

printf ( "subfam. .. invalid format specified for freqzero2\n" ) ; 
exit () ; 

} /* end switch */ 

/* */ 

/* MULTIPLICATION */ 

/* */ 

/* create space for the output array */ 

/* */ 

s3= (COMPLEX* ) calloc (n, sizeof (COMPLEX) ) ; 
if (s3==NULL) { 

printf ( "subfam. .. insufficient space to allocate s3\n"); 
exit ( ) ; 

} 

if ( (conjl [0] ==' y' ) && (conj2 [0] ==' y' ) ) { /* multiply conjugate si times conjugate s2 

for (i=0; i<n; i++) { 

s 3 [ i ] . r=sl [i] . r*s2 [i] . r-sl [i] .i*s2[i] . i; 
s 3 [ i ] . i= (-1) * (si [i] . i*s2 [i] . i+sl [i] .i*s2 [i] .r) ; 

} 

} 

if ( (con j 1 [0 ] ==' y' ) && (con j2 [0 ] ==' n' ) ) { /* multiply conjugate si times s2 */ 

for (i=0; i<n; i++) { 

s 3 [ i ] .r-sl [i] . r*s2 [i] . r+sl [i] . i*s2 [i] . i; 
s3 [i] . i=sl [ i] . r*s2 [i] . i-sl [i] . i*s2 [i] . r; 

} 

} 

if ( (conjl [0] = =/ n' )&& (conj2 [0]==' y' ) ) { /* multiply si times conjugate s2 */ 

for (i=0; i<n; i++) { 

s 3 [ i ] .r=sl [i] . r*s2 [i] .r+sl [i] . i*s2 [i] . i; 
s3[i] .i=sl[i] .i*s2[i] .r-sl[i] .r*s2[i] .i; 

} 

} 

if ( (conjl [0]= =/ n' )&& (conj2 [0]= =/ n' ) ) { /* multiply si times s2 */ 
for (i=0; i<n; i++) { 

s 3 [ i ] . r=sl [i] . r*s2 [i] .r-sl [l] .i*s2[i] . i; 
s3 [i] . i=sl [i] . r*s2 [i] .i+sl [i] .i*s2[i] . r; 
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} 

\ 

/* take an FFT of n-point s3 */ 

/* set values in arguments sent to fft */ 

direct ion=l; / * request a forward transform time-to-freq * 

norm=0 ; 

fft (s3, n, direction, norm) ; /* results returned in s3 */ 

/* */ 

/* OUTPUT */ 

/* */ 

/* take magnitude and phase of s3 and place in the output file * 
/* */ 

switch (ioflag[2]) { 

case 'a' : /* outputfile is an ascii file */ 

if (( ofp = fopen (out file, "w" ) ) == NULL ) { 

print f ("subf am. . .unable to open ascii output file\n" ) ; 
exit () ; 

} 

/* write out data to the ascii file */ 

for (i=0; i<n; i++) { 

outputmag=sqrt ( (s3 [ i ] . r*s3 [ i] . r ) + (s3 [i] . i*s3 [ i] . i) ) ; 
outputphase=atan (s3 [i ] . i/s3 [ i] . r) ; 

fpnntf (ofp, "%-16. 6e %-16. 6e\n", outputmag, outputphase) ; 

} 

/* close the output file */ 

f close (ofp) ; 

break; 

case ' b' : /* outputfile is a binary file */ 

if (( ofp = fopen (outfile, "wb" ) ) == NULL ) { 

printf (" subf am. . .unable to open binary output file\n" ) ; 
exit ( ) ; 

} 

/* write out data to the binary file */ 

for (i=0; i<n; i++) { 

outputmag=sqrt ( (s3 [i] . r*s3 [i] . r) + (s3 [i] . i*s3 [i] .i)); 
outputphase=atan (s3 [i] . i/s3 [i] . r) ; 

*y=outputmag; 

if (f write (y, sizeof (*y) , 1, ofp) ! =1) { 

printf (’’subfam. . .error while writing output file\n" ) ; 
fclose (ofp) ; 
exit () ; 

} 

*y=outputphase ; 

if ( f write (y, sizeof ( *y) , 1, ofp) 1=1) { 

print f ( "subfam. .. error while writing output file\n" ) ; 
fclose (ofp) ; 
exit () ; 

} 

} /* end if i=0 ... */ 

fclose (ofp) ; 

break; 

default: /* invalid format specified for outputfile */ 

print f ("subfam. .. invalid format specified for outputf ile\n" ) 
exit () ; 

} /* end switch ioflag[2] */ 
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LIST OF ABBREVIATIONS 



BPSK 

FAM 

FFT 

FSM 

PAM 

QAM 

QPSK 

SCF 

SSCA 

SSPI 

SUB FAM 


binary phase shift keying 
FFT accumulation method 
fast fourier transform 
frequency smoothing method 
pulse amplitude modulation 
quadrature amplitude modulation 
quadrature phase shift keying 
spectral correlation function 
strip spectral correlation algorithm 
Statistical Signal Processing, Inc. 
sub -FFT accumulation method 
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